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ABSTRACT 


This  report  details  the  development  and  evaluation  of  a  new  direction  finding 
estimator,  the  Spread  Maximum  Likelihood  algorithm,  targeted  for  high  latitude  HF  prop¬ 
agation  conditions.  Of  particular  concern  are  patches  of  enhanced  electron  density  and 
associated  instabilities  in  the  ionosphere  which  drift  across  the  polar  cap  during  dark¬ 
ness  causing  scattering  of  the  propagated  signal  over  a  range  of  azimuth  and  elevation 
directions  (the  signals  are  spatially  spread).  The  new  estimator  incorporates  a  spatial 
spreading  model  allowing  it  to  simultaneously  deal  with  both  patch  scattered  signals  and 
signals  propagated  by  more  normal  propagation  mechanisms,  as  well  as  to  distinguish  be¬ 
tween  the  two  types.  Evaluation  of  the  new  estimator  using  both  simulation  and  collected 
data  show  it  to  be  considerably  superior  to  both  conventional  and  modern  superresolution 
approaches  for  high  latitude  propagation  conditions. 

RESUME 

Ce  rapport  decrit  une  technique  pour  l’estimation  de  Tangle  de  gisement  des 
ondes  HF  ainsi  que  son  evaluation.  L’algorithme  intitule  “vraisemblence  de  l’etalement 
maximal”  s’applique  aux  ondes  HF  transmisse  dans  les  latitudes  elevee.  II  s’applique 
particulierement  aux  zones  de  l’ionosphere  ou  des  perturbations  associees  4  une  densite 
d ’electrons  plus  elevee  se  produisent.  Ces  zones  derivent  au-dessus  de  la  calotte  polaire 
pendant  la  nuit,  ce  qui  entraine  une  diffusion  du  signal  propage  sur  une  gamme  d’azimuts 
et  de  directions  d’elevation,  ce  qui  veut  dire  que  les  signaux  subissent  un  etalement  spatial. 
La  technique  d’estimation  comprend  un  module  d’etalement  spatial  qui  permet  de  traiter 
les  signaux  se  propageant  a  traver  les  zones  instables  ainsi  que  les  signaux  propages  par 
des  mecanismes  plus  conventionnels  et  de  faire  la  distinction  entre  les  deux.  L’evaluation 
menee  avec  des  donnees  fictives  et  reelles  a  demontre  que  la  methode  proposee  est  superieur 
aux  approches  conventionnelles  et  modernes  de  superresolution  pour  la  propagation  en 
haute  latitude. 
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EXECUTIVE  SUMMARY 


The  requirement  exists  to  improve  the  accuracy  of  strategic  high  latitude  HF 
direction  finding  (DF)  systems.  In  the  past,  poor  DF  accuracy  derived  from  Arctic  mea¬ 
surements  has  led  to  low  confidence  in  high  latitude  sites  despite  the  strategic  relevance 
of  these  sites  for  transmitter  geolocation. 

The  biggest  problem  appears  to  be  patches  of  enhanced  electron  density  and 
associated  instabilities  in  the  F  layer  of  the  ionosphere  which  drift  across  the  polar  cap 
during  darkness.  These  patches  scatter  the  signal  over  a  range  of  azimuth  and  elevation 
directions  (the  signals  are  spatially  spread)  which  are  very  different  from  the  expected 
bearing  given  the  transmitter’s  location. 

Although  the  bearings  of  these  spatially  spread  signals  provide  no  indication  of 
the  true  bearing  of  the  transmitter,  there  is  evidence  that  sporadic-E  propagation  (i.e. 
reflections  from  the  E  layer  of  the  ionosphere)  may  sometimes  be  supported  during  these 
conditions  and  that  this  propagation  mode  yields  the  desired  signal  bearing.  Unfortu¬ 
nately,  conventional  DF  approaches  estimate  the  bearing  of  the  most  dominant  propaga¬ 
tion  mode  which  will  generally  be  the  patch  scattered  signals.  Newer  superresolution  DF 
algorithms,  which  can  simultaneously  DF  signals  from  several  directions,  tend  to  become 
overloaded  by  the  range  of  directions  occupied  by  the  scattered  signals.  Hence  there  is 
a  need  for  a  new  approach  which  is  capable  of  dealing  with  this  unique  high-latitude 
phenomena. 

In  this  report,  a  new  high  latitude  HF  DF  algorithm  is  introduced.  The  theo¬ 
retical  development  of  this  new  algorithm  is  based  on  maximum  likelihood  principles  and 
the  key  feature  is  that  a  model  of  the  spatial  spreading  is  incorporated.  The  performance 
of  this  new  algorithm  -  the  spread  maximum  likelihood  (SML)  method  -  is  evaluated 
using  both  simulated  data  and  off-air  data  collected  at  the  Arctic  site  CFS  Alert.  This 
evaluation  included  the  effects  of  noise,  the  amount  of  spatial  spreading,  and  the  shape  of 
the  spread  region.  It  also  included  an  evaluation  of  the  ability  to  detect  a  weaker  point- 
source  (no  spatial  spreading)  signal  in  the  presence  of  stronger  spatially  spread  signals. 
The  results  of  this  evaluation  show  that  the  new  algorithm  is  not  only  capable  of  dealing 
with  spatially  spread  signals  but  also  provides  information  on  the  amount  of  spreading. 
This  extra  information  is  invaluable  since  it  provides  a  useful  means  of  qualifying  signal 
accuracy.  For  example,  a  large  amount  of  spatial  spreading  indicates  a  patch  scattered 
signal  which  can  be  ignored,  while  little  spreading  indicates  the  signal  was  propagated 
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by  a  more  conventional  propagation  mode  and  more  likely  yields  the  correct  transmitter 
bearing. 

As  is  already  evident,  the  main  advantages  of  the  SML  algorithm  are  that,  for 
high  latitude  HF  conditions,  it  provides  more  useful  information  about  the  received  sig¬ 
nals,  and  it  is  considerably  more  accurate  than  both  conventional  and  modern  superreso¬ 
lution  DF  approaches  when  propagation  conditions  are  poor  (i.e.,  when  patch  scattering 
occurs).  For  more  benign  conditions,  the  performance  of  the  SML  algorithm  is  at  least  as 
good  as  the  other  approaches.  For  these  reasons,  serious  consideration  should  be  given  to 
its  implementation  at  high  latitude  DF  sites.  * 

The  main  disadvantage  of  the  SML  algorithm  is  that  it  is  computationally  very 
intensive  which  makes  it  an  order  of  magnitude  or  more  slower  than  other  approaches.  For 
operational  purposes,  this  problem  could  be  overcome  by  using  the  SML  algorithm  only 
when  patch  scattering  conditions  are  suspected,  and  using  faster  DF  methods  otherwise. 

Future  research  should  focus  on  developing  realtime  implementations  of  the  SML 
algorithm.  This  is  a  possibility  which  could  be  realizable  within  the  next  few  years  given 
the  ever  increasing  speed  of  modern  computing  power. 
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1.0  INTRODUCTION 


The  requirement  exists  to  improve  the  accuracy  of  strategic  HF  direction  finding  sys¬ 
tems,  particularly  in  the  Arctic.  In  the  past,  poor  DF  accuracy  derived  from  Arctic 
measurements  has  led  to  low  confidence  in  high  latitude  sites  despite  the  strategic  rele¬ 
vance  of  these  sites  for  transmitter  geolocation. 

The  biggest  problem  appears  to  be  patches  of  enhanced  electron  density  and  associated 
instabilities  in  the  F  layer  of  the  ionosphere  which  drift  across  the  polar  cap  during 
darkness  in  a  roughly  antisunward  direction  at  speeds  ranging  from  a  few  hundreds  to 
over  one  thousand  meters  per  second  [1].  These  patches  can  cause  scattering  of  a  signal 
from  azimuth  directions  which  are  very  different  from  the  true  bearing  of  the  transmitter. 

One  avenue  of  investigation  being  pursued  is  the  development  of  new  DF  algorithms 
which  are  better  matched  to  the  high  latitude  HF  signal  environment.  Currently  popular 
DF  algorithms  assume  that  the  incoming  signal  can  be  modeled  as  a  point  source,  or 
equivalently,  the  incoming  signal  has  a  planar  wavefront.  This  is  reasonable  if  the  size  of 
the  transmission  source  is  extremely  small  relative  to  its  range;  the  size  of  the  DF  array  is 
also  small  relative  to  the  transmitter  range;  the  ionosphere  acts  as  a  perfect  or  near-perfect 
reflector;  and  local  site  multipath  can  be  ignored.  Unfortunately  at  high  latitudes,  during 
periods  where  scattering  from  large  moving  patches  is  occurring,  the  received  signals  are 
distributed  over  a  range  of  bearings  in  both  azimuth  and  elevation.  Current  DF  algorithms 
handle  this  phenomena  by  splitting  it  up  into  several  signals  coming  from  the  same  area 
of  the  ionosphere. 

Although  representing  a  signal  scattered  from  a  region  of  the  ionosphere  (called  a 
spread-source  signal  in  this  report)  by  several  point-source  signals  is  still  informative,  it 
degrades  the  ability  of  the  algorithm  to  detect  weaker  signals  that  may  also  be  present.  For 
example,  previous  studies  have  indicated  that  the  E  layer  of  the  ionosphere  acts  as  a  good 
reflector  when  it  supports  signal  propagation  (i.e.  sporadic-E)  between  the  transmitter 
and  receiver.  Hence,  it  is  important  that  a  DF  algorithm  be  able  to  detect  and  determine 
the  bearing  of  a  weaker  sporadic-E  reflected  signal  in  the  presence  of  one  or  more  stronger 
spread  signals. 

In  attempts  to  meet  these  requirements,  a  number  of  DF  algorithms  were  developed 
at  McMaster  University  [2].  Although  these  algorithms  showed  promise,  in  their  develop¬ 
ment  a  strong  effort  was  made  to  keep  the  processing  considerations  reasonable  and  this, 


unfortunately,  resulted  in  algorithms  which  were  only  able  to  deal  with  a  single  signal. 

To  overcome  the  single  signal  limitation,  processing  considerations  were  abandoned  in 
favour  of  developing  the  most  accurate  approach  possible  capable  of  dealing  with  mul¬ 
tiple  spread  signals.  To  this  end,  the  maximum  likelihood  approach  was  chosen,  since 
this  approach  leads  to  optimal  estimators  when  all  known  information  about  the  signal 
environment  is  taken  into  account.  The  result  is  a  new  estimator  called  the  Spread  Max¬ 
imum  Likelihood  Estimator  whose  development  and  testing  is  the  subject  of  the  rest  of 
this  report. 

The  layout  of  the  rest  of  the  report  is  as  follows.  In  Section  2,  the  maximum  likelihood 
method  is  introduced  followed  by  the  development  of  the  specific  signal  model  used  by 
the  spread  maximum  likelihood  (SML)  variant.  In  Sections  3-5,  a  detailed  description  of 
the  SML  algorithm  is  provided,  including  the  initial  value  assessment  in  Section  3,  the 
refinement  procedure  in  Section  4,  and  the  step-by-step  algorithm  summary  in  Section 
5.  In  Section  6,  the  optimum  settings  for  various  algorithm  parameters  is  investigated 
through  simulation.  Comparisons  are  also  made  with  the  MUSIC  DF  algorithm  (an 
algorithm  that  is  commonly  used  in  advanced  DF  processing)  using  simulated  data  and 
off-air  data  collected  with  a  12-channel  DF  system  in  the  Arctic.  Finally,  in  Section  7  the 
conclusions  and  recommendations  are  presented. 
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2.0  SPREAD  MAXIMUM  LIKELIHOOD  METHOD  OVERVIEW 


2.1  Maximum  Likelihood  Estimation 


A  successful,  albeit  often  computationally  intensive,  approach  to  estimation  is  based 
on  the  maximum  likelihood  method.  Essentially  the  idea  is  to  find  the  most  likely  state  of 
a  signal  process  given  a  set  of  measurement  observations  made  of  the  process.  Since  this 
is  a  statistical  approach,  the  method  applies  to  cases  where  either  the  signal  generation 
is  a  random  process,  or  the  measurements  have  been  corrupted  by  additive  noise  (which 
is  a  random  process),  or  both. 


Assuming  the  random  processes  are  all  normally  Gaussian  distributed,  and  measure¬ 
ments  are  made  using  N  sensors,  then  the  associated  probability  density  function  is  given 
by  [3] 


/(x  0,Xi, 


x„  A  - _ i _  -traceax-M^c-^x-M)) 

K~l)  KdetC]*  6 


(1) 


where  the  superscript  H  denotes  the  conjugate-transpose  operation,  the  vectors  x0, ...,  *k-i 
represents  the  random  complex  measurement  data  associated  with  all  N  sensors  for  time 
instances  t  —  t0,ti,  1  as  defined  by 


xfc  = 


x0(k) 

Xi(k) 

x2(k) 


for  0  <  k  <  K. 


xN-i(k) 


(2) 


Additionally,  the  matrix  X  is  merely  a  compact  way  of  representing  all  K  measurement 
vectors  and  is  given  by 

X  =  [x0,xi,  (3) 

the  matrix  M  represents  the  corresponding  mean  values  of  the  measurements,  and  C 
is  the  N  y.  N  covariance  matrix  describing  the  correlations  among  sensors.  The  exact 
definitions  of  the  matrices  M  and  C  depend  on  how  the  above  density  function  is  applied 
to  the  particular  estimation  problem  to  be  solved,  as  will  be  seen.  Once  this  equation 
has  been  properly  set  up,  the  maximum  likelihood  solution  is  then  found  by  maximizing 
the  probability  density  function  by  making  the  appropriate  choice  of  M  and/or  C  from 
among  the  allowable  choices.  The  optimum  choices  of  M  and/or  C  can  then  be  related 
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back  to  the  signal  parameters  of  interest. 

For  direction  finding  purposes  Xi, xk-i  are  taken  to  represent  the  complex  baseband 
outputs  from  an  array  of  N  antennas.  The  definitions  for  M  and  C  depends  on  the 
assumptions  made  about  the  signals.  The  two  most  commonly  used  assumptions  are 
either 

1.  the  signal  is  deterministic  but  unknown,  or, 

2.  the  signal  is  stochastic  with  a  zero-mean  Gaussian  distribution. 

Choosing  the  first  assumption,  the  data  can  be  decomposed  as 

Xn(k)  =  mn(k)  +T}n(k)  (4) 

where  mn(k )  is  the  exact  (but  unknown)  signal  value  and  rjn(k)  is  the  value  of  the  noise 
corrupting  the  measurement.  From  this  expression  it  follows  that  the  elements  of  M  are 
given  by  mn(k),  and  the  elements  of  C  then  represent  the  covariance  of  the  noise  processes 
only,  as  defined  by 

Cmn  =  E[qm(k)r)n(ky]  for  0  <  k  <  K  and  0  <  m,n  <  N.  (5) 

The  approach  based  on  the  first  assumption  is  called  the  deterministic  maximum  likelihood 
method. 

Choosing  the  second  signal  assumption  leads  to 

mn(k )  =  0  (6) 

and 

Cmn  =  E[xm(k)xn(ky]  for  0  <  k  <  K  and  0  <m,n<N  (7) 

since  the  sum  of  two  or  more  zero-mean  Gaussian  processes  is  still  a  zero-mean  Gaus¬ 
sian  process.  The  approach  based  on  the  stochastic  assumption  is  called  the  stochastic 
maximum  likelihood  method. 

Comparing  these  two  approaches  reveals  that  the  main  differences  are  the  constraints 
placed  on  the  estimated  complex  signal  amplitude(s).  The  deterministic  signal  assumption 
(usually)  leads  to  no  constraint  while  the  stochastic  signal  assumption  leads  to  a  Gaussian 
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distribution.  Given  that  man-made  signals  are  bounded  (finite  power)  and  often  have  a 
Gaussian-like  distribution,  the  second  assumption  generally  results  in  better  estimator 
performance  [4]  (see  also  [5]  for  a  comparison  of  the  deterministic  approach  with  better 
modeled  methods).  An  exception  to  this  might  be  the  case  of  constant  modulus  signals 
(e.g.,  fm,  psk,  fsk,  etc.)  since  a  constant  amplitude  constraint  could  more  easily  be 
added  to  the  deterministic  approach  than  the  stochastic  approach.  However,  for  the  HF 
signals  considered  in  this  report,  the  dynamic  nature  of  the  ionosphere  tends  to  destroy 
the  constant  modulus  properties  of  the  received  signal,  so  that  the  stochastic  maximum 
likelihood  approach  was  followed  (in  fact  the  deterministic  approach  was  tried  but  was 
not  found  to  be  nearly  as  successful). 


2.2  The  Cost  Function 

In  this  section,  the  stochastic  maximum  likelihood  method  is  further  refined.  The 
central  key  to  this  refinement  is  the  development  of  a  signal  plus  noise  model  which 
accurately  reflects  the  high  latitude  HF  environment.  This  involves  accounting  for  signals 
which  are  diffusely  reflected  or  scattered  off  a  highly  disturbed  region  of  the  ionosphere 
with  an  angular  extent  of  several  degrees  or  more  in  both  azimuth  and  elevation,  hence 
the  name  Spread  Maximum,  Likelihood  Estimator.  By  comparison,  the  standard  stochastic 
maximum  likelihood  approach  considers  only  point-sources  signals. 

Before  beginning  the  model  development,  the  probability  density  function  (1)  can  be 
simplified  using  (6)  to  give 

/(x„,x1,...,xK_1)  =  e-trace(x*c-'x)  (8) 

by  using  the  stochastic  signal  model.  The  objective  remains  the  same  as  before,  that  is, 
find  the  unknown  covariance  matrix  C  which  maximizes  this  expression.  This  is  equivalent 
to  maximizing  the  cost  function  Lc  =  ln/(x0,Xi,  ...,x.K_i)  or, 

Lc  =  —NK  ln(7r)  -  K  ln(det  C)  -  trace(X* C_1X) .  (9) 

Since  the  addition  or  multiplication  by  a  constant  value  has  no  effect  on  the  maximization, 
the  cost  function  can  be  simplified  to 

L  =  —  ln(det  C)  -  trace(RC_1)  (10) 
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where  R  is  the  data  covariance  matrix  constructed  from  X  using 

R=ixx».  (11) 

2.3  Modeling  the  Signal  Environment 

The  modeling  aspect  comes  into  play  when  C  is  being  selected.  A  model  is  used  to 
generate  C  based  on  input  modeling  parameters  such  as,  for  example,  the  number  of 
signals,  signal  bearings,  signal  amplitudes/powers,  and  noise  powers.  For  this  reason, 
C  is  referred  to  as  the  model  covariance  matrix  in  this  report.  Once  the  model  has 
been  setup,  a  common  approach,  and  the  approach  employed  here,  is  to  choose  some 
initial  model  parameters,  generate  C,  and  then  determine  the  cost  function  value  L.  The 
model  parameters  are  then  successively  refined  and  C  recomputed  until  the  maximum 
cost  function  value  has  been  achieved.  The  model  parameter  values  corresponding  to  this 
maximum  value  are  then  taken  to  be  the  optimum  or  maximum  likelihood  estimates. 

The  particular  model  used  to  generate  the  model  covariance  matrix  depends  on  many 
factors  including  the  transmitter (s)  and  receiver  characteristics,  the  signal  propagation 
environment,  the  noise  sources,  and  so  on.  One  way  to  set  this  model  up,  is  to  consider  the 
generation  of  synthetic  data  which  imitates  the  collected  data,  and  then  use  this  synthetic 
data  to  determine  C  according  to 

C  =  ^YYh  (12) 

where  Y  is  the  matrix  of  model  data  and  has  the  same  form  and  dimensions  as  X. 

Proceeding  in  this  manner,  it  becomes  necessary  to  make  certain  assumptions  about 
the  environment  in  order  to  simplify  the  derivations.  Initially,  the  following  ideal  assump¬ 
tions  can  be  made: 

1.  M  incoming  point-source  signals  arrive  from  M  different  directions. 

2.  The  noise  is  additive  and  Gaussian. 

3.  The  receiver  channels  are  linear  and  perfectly  matched. 

4.  The  receive  antennas  are  isotropic  and  located  at  the  same  height. 

5.  There  is  no  coupling  among  the  receive  antennas  or  with  the  local  environment. 
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Modifications  to  these  assumptions  will  then  be  introduced  as  appropriate. 

Based  on  these  simplifying  assumptions,  the  received  signal  can  be  decomposed  as 

Y  =  Ya+Y2  +  ...+YM  +  N  (13) 

where  the  matrices  Yi,  Y2, ...,  YM  represent  the  model  data  for  the  M  individual  signals, 
and  N  is  the  modeled  noise.  The  matrix  Ym,  for  0  <  m  <  M,  can  be  defined  in  vector 
form  as 

Ym  =  ema£  (14) 

where  em  is  the  array  response  (or  steering)  vector  for  the  mth  signal,  and  am  is  the 
corresponding  signal  amplitude  vector.  The  definitions  for  the  elements  of  the  steering 
vector  are  given  by 


1 

Vn 


eix(Xo  sin  C0S  ^m+yo  COS  4>m  COS  tpm) 

ej  ^  (rri  sin  <£m  cos  ^m+yi  cos  <j>m  cos  ipm ) 


sin</>m  cos  ipm ~\~Vn — l  cos  </>m  cos  ipm) 


(15) 


where  A  is  the  signal  wavelength,  xn  and  yn  are  the  Cartesian  coordinates  for  antenna 
n  (with  the  phase  center  of  the  array  located  at  the  origin),  <f)m  is  the  azimuth  angle  of 
the  mth  signal  measured  clockwise  with  respect  to  the  Y-axis  of  the  coordinate  system, 
and  'ipm  is  the  elevation  angle  measured  with  respect  to  the  X-Y  plane  (the  ground).  The 
definition  for  the  elements  of  the  signal  amplitude  vector  is  given  by 


flm(O)* 

Om(l)‘ 

am(K  -  1)* 


(16) 


where  am(0),am(l),  ...,am(K  —  1)  are  the  received  complex  amplitudes  of  the  mth  signal 
for  time  instances  to,ti, ...,  t^-i- 

Analytically,  the  preceding  approach  to  generating  the  model  covariance  matrix  is 
useful,  however,  it  requires  too  many  model  parameters.  For  example,  the  generation  of 
Y  requires  selecting  the  azimuth  angle,  elevation  angle,  and  complex  amplitudes  for  each 
signal  plus  the  corresponding  complex  noise  values  for  each  sensor.  With  the  exception 
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of  the  angles,  this  must  be  repeated  for  each  time  instance  resulting  in  a  total  number 
of  2 M  +  2  (M  +  N)K  real  model  parameters  (remembering  that  one  complex  parameter 
is  equivalent  to  two  real  parameters).  Additionally  the  noise  and  signal  amplitudes  must 
also  be  chosen  so  that  statistically  they  have  complex  Gaussian  distributions. 

A  more  practical  approach  is  to  rewrite  the  model  covariance  as  a  sum  of  the  noise 
covariances,  the  signal  covariances,  and  the  signal  cross-covariances,  or 

M  M—\  M 

C  =  <r2C„+£Cmm+£  £  (Qmn  +  Qm«)  (17) 

771=1  771=1  71=771+1 

where  a2  is  the  noise  power,  Cv  is  the  normalized  noise  covariance  matrix  (trace  C,,  =  1), 
Cmm  is  the  signal  covariance  matrix  for  the  mth  signal,  and  Qmn  is  the  signal  cross¬ 
covariance  matrix.  The  generation  of  these  matrices  is  discussed  in  the  following  para¬ 
graphs. 

The  noise  covariance  matrix  Cv  is  assumed  to  be  known  a  priori  and  will  not  be 
considered  as  part  of  the  estimation  process.  The  determination  of  Cv  can  be  done  either 
through  theoretical  statistical  considerations,  or  through  measurements.  For  example,  if 
the  noise  is  known  to  be  white  Gaussian  in  nature  with  equal  but  uncorrelated  amplitudes 
in  each  channel,  then  E[r}m{k)r)m(k)*]  =  a2  and  E[rjm(k)r]n(k)*]  =  0  for  0  <  m,n  <  N 
and  m^n,  hence 

c,  =  (18) 

where  Ijv  is  the  N  x  N  identity  matrix.  Alternatively,  if  data  measurements  can  be  taken 
when  no  signals  are  present,  then 

^  XXff 

v  ~  trace  (XXff) '  ^ 

More  elaborate  procedures  could  be  developed  involving  a  number  of  measurement  sets 
with  the  same  noise  environment  but  different  signal  directions,  however  the  development 
of  this  kind  of  approach  is  beyond  the  scope  of  this  report.  It  suffices  to  say  that  joint 
estimation  of  both  the  noise  and  signal  characteristics  should  be  avoided  if  possible  since 
it  leads  to  poorer  accuracy. 


The  signal  covariance  matrix  Cmm  can  be  defined  in  terms  of  the  model  data  as 


C—  __v 
mm  —  x  m  x  m 
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which  can  be  further  expanded  in  terms  of  the  component  amplitude  and  steering  vectors 
as  R 

p  _  amam  H  /'oi\ 


In  a  similar  fashion,  the  cross-covariance  matrix  Qmn  can  also  be  defined  and  expanded 
to  get 


Q 


ran 


H 

n 


K 


6me 


H 

71 


for  n  >  m. 


(22) 


At  this  point,  the  number  of  model  parameters  has  not  been  changed  (except  for  the 
noise  covariances  which  are  now  assumed  to  be  known).  Simplifications  can  be  introduced 
to  the  expressions  for  Cmm  and  Qmn  by  recognizing  that  the  term  a^. /K  is  the  average 
signal  power  over  the  given  time  interval  and  can  be  replaced  by  the  single  parameter  sm. 
Equation  (21)  can  then  be  rewritten  as 


C 


mm 


Sm^m^rn- 


(23) 


Additionally,  defining  the  complex  correlation  coefficient 


=  am^n 

\&m  |  |*^n| 


where  |a|  =  (aHa)z,  then  Qmn  can  be  rewritten  as 


(24) 


Qmn  —  Pmn(.SmSn)  2  ' 


(25) 


Hence  the  K  complex  amplitude  model  parameters  for  each  signal  have  been  replaced  by 
a  single  power  parameter  plus  up  to  M  —  1  complex  correlation  coefficients.  Also  counting 
the  bearing  angles  required  for  each  signal  (4>m  and  ipm),  and  the  noise  power,  then  the 
number  of  real  model  parameters  is  3 M  +  M(M  - 1)  + 1  which  is  a  considerable  reduction 
in  parameters  when  K  M  (which  will  normally  be  true  in  practice). 

A  further  major  reduction  in  the  number  of  parameters  is  possible  if  the  correlations 
pmn  are  known.  For  the  HF  skywave  environment,  signals  from  different  transmitters  will 
be  uncorrelated.  Signals  originating  from  the  same  transmitter  but  reflected  off  different 
layers  of  the  ionosphere  will  also  be  uncorrelated  due  to  the  large  path  length  differences 
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usually  encountered  (i.e.,  the  path  delay  time  differences  will  be  greater  than  the  inverse 
bandwidth  of  the  signal).  Hence  one  can  assume  that 

Pmn  =  0  for  0  <  m  <  n  <  M  (26) 

which  eliminates  the  third  term  in  (17)  and  reduces  the  number  of  real  model  parameters 
to  3  M  +  1. 

A  violation  of  the  assumption  of  uncorrelated  signals  is  the  case  of  local  site  multipath. 
However,  like  the  problem  of  determining  the  noise  correlations,  including  the  signal 
correlations  in  the  estimation  process  is  highly  undesirable  (although  in  the  standard 
stochastic  ML  approach  this  is  done  [6]).  Hence  it  is  assumed  that  either  the  receiver  site 
is  well  chosen,  or  signal  correlations  can  be  determined  independently  and  corrected  as 
done,  for  example,  in  [7]. 

Up  until  this  point,  only  incoming  point-source  signals  have  been  considered.  For  HF, 
this  assumes  perfect  reflections  off  the  ionosphere,  which  is  a  reasonable  assumption  for 
benign  atmospheric  conditions.  At  high  latitudes,  however,  the  ionosphere  can  sometimes 
become  highly  disturbed  [8]  resulting  in  behaviour  which  is  not  consistent  with  an  ideal 
reflector.  For  example,  measurements  have  shown  that  during  these  disturbed  conditions 
Doppler  shifts  of  up  to  45  Hz  or  more  can  occur  in  conjunction  with  large  eigenvalue 
spreading  in  the  measured  data  covariance  matrix  [9].  These  measurements  suggest  that 
the  transmitter  signal  is  scattered  from  highly  dynamic  regions  of  the  ionosphere  yielding 
the  large  Doppler  spreads.  These  measurements  also  suggest  that  the  regions  are  suffi¬ 
ciently  wide  (angular  spreads  of  several  degrees  or  more)  that  signal  decorrelation  occurs 
across  the  region  causing  the  large  eigenvalue  spreading. 

This  interpretation  is  also  supported  by  HF  propagation  research  work  carried  out 
at  CRC  [10].  This  work  suggests  that  the  most  likely  mechanism  for  the  ionospheric 
disturbances  is  moving  patches  of  enhanced  electron  density  with  sizes  greater  than  several 
kilometers  which  drift  across  the  polar  cap  during  darkness.  As  these  patches  progress, 
they  give  rise  to  short-lived  small-scale  electron  irregularities  on  the  scale  of  less  than  one 
kilometer.  These  irregularities  become  very  elongated  along  the  earth’s  magnetic  field 
lines  and  tend  to  scatter  incident  radio  waves  through  360°  in  azimuth  with  little  or  no 
scattering  in  elevation.  Given  the  random  and  short-lived  nature  of  these  irregularities, 
the  patch  acts  as  a  temporal  and  spatial  incoherent  radio  scattering  region. 

The  interpretation  of  conditions  encountered  at  high  latitudes  leads  to  a  modified 
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Figure  1:  Grid  of  point-source  signals  used  to  represent  a  spread  signal.  Each  grid  node 
represents  location  and  power  of  each  corresponding  source  signal. 


description  for  the  signal  covariance  matrix  Cmm.  To  begin  with,  the  scattering  region 
is  divided  into  a  grid  where  each  section  in  the  grid  is  made  small  enough  that  it  can  be 
accurately  represented  by  a  single  perfectly  reflected  point-source  signal.  Figure  1  shows 
an  example  plotting  the  power  of  the  point-source  signals  versus  azimuth  and  elevation 
angle  for  a  single  region.  The  power  towards  the  edges  of  the  modeled  region  is  tapered 
according  to  some  predetermined  taper  function  which  is  chosen  to  best  represent  the 
actual  environment.  Given  that  signals  scattered  from  different  parts  of  the  same  region 
are  uncorrelated,  the  model  covariance  matrix  can  be  generated  according  to  the  sum  of 
the  individual  covariances  of  the  component  point-source  signals,  or 

Mg 

^ mm  y  l  (27) 

i=l 

where  Mg  is  the  number  of  component  point-source  signals  in  the  grid,  emi  represents 
the  steering  vectors  of  the  grid  signals,  /?m($,  '3')  is  the  taper  function  which  is  positive 
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real-valued  and  normalized  so  that 

Mg 

I =  1  (28) 

2—1 

and  $mi-  and  are  the  normalized  grid  angles.  The  normalized  grid  angles  are  defined 
by 

(29) 

(30) 

where  (j>m  and  are  the  respective  azimuth  and  elevation  bearing  angles  of  the  center 
of  the  grid,  <j)mi  and  are  the  respective  azimuth  and  elevation  angles  of  the  compo¬ 
nent  point-source  signals,  and  A^m  and  A^m  are  the  signal  spread  parameters  defining 
the  angular  extent  of  the  spread  signal  (width  of  the  grid)  in  azimuth  and  elevation, 
respectively. 

In  order  to  be  mathematically  usable,  the  taper  function  must  vary  smoothly  as  a 
function  of  the  azimuth  and  elevation  spread  parameters,  and  the  differentiated  taper 
function  must  be  defined  for  all  valid  values  of  the  spread  parameters.  For  example,  in 
a  rectangular  grid  with  uniform  fixed  spacings  and  equal  powers  for  all  individual  point- 
source  signals,  the  taper  function  changes  in  discrete  steps  as  each  new  column/row  of 
point-source  signals  is  introduced  or  removed  when  the  azimuth/elevation  spread  changes 
sufficiently.  The  choice  of  taper  function  and  grid  layout  (since  this  affects  the  taper 
function)  is  discussed  in  Section  4.4. 

Using  (17),  (27)  and  incorporating  the  simplifications  to  the  model  covariance  gener¬ 
ation  already  described,  then 

M  Mg 

C  =  O  c„  +  Sm  ^ mij^mi^rni  (31) 

m=  1  i=l 

where  the  noise  power  <r2,  signal  power  sm,  azimuth  bearing  <f)m,  elevation  bearing  ipm, 
azimuth  spread  A^m,  and  elevation  spread  A^m  are  the  only  input  values  which  are 
explicitly  adjusted  in  order  to  maximize  the  cost  function  L. 


^  $mi 


Vlr  •  — 

^  mr  — 


i<Pm 

An  -  A 
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3.0  MODEL  PARAMETER  ESTIMATION 


The  previous  section  examined  the  generation  of  the  model  covariance  matrix  and  its 
evaluation  via  the  cost  function.  In  this  section,  an  overview  of  the  model  parameter 
estimation  procedure  is  provided  along  with  details  on  determining  the  initial  model 
parameter  estimates.  Details  on  the  parameter  refinement  procedure  are  provided  later 
in  Section  4.0. 

A  major  problem  with  determining  the  optimum  model  parameters  in  maximum  likeli¬ 
hood  problems  is  that  it  is  usually  difficult  or  impossible  to  develop  a  closed-form  solution. 
The  standard  approach  is  to  employ  some  form  of  search  procedure  to  find  the  solution. 
An  exhaustive  search  of  all  possible  model  parameter  values  could  be  performed  to  do  this, 
however,  with  more  than  a  very  small  number  of  model  parameters,  this  is  a  prohibitively 
time-consuming  process.  Faster  approaches  such  as  gradient  ascent/descent  techniques 
are  useful,  but  require  the  initial  model  parameter  to  be  reasonably  close  to  the  optimum 
solution,  otherwise  the  search  may  converge  on  a  false  solution.  An  approach  which  has 
been  used  successfully  in  the  past,  and  works  reasonably  well  with  uncorrelated  signals,  is 
the  Alternating  Projection  Maximization  (APM)  method  [11].  Although  heavily  modified, 
the  approach  employed  here  is  conceptually  similar. 

The  basic  procedure  is  to  build  up  the  model  covariance  matrix  one  signal  at  a  time. 
In  the  very  first  stage  an  estimate  of  the  noise  power  is  made  assuming  the  data  consists 
of  noise  only  (Section  3.1).  Using  this  as  the  starting  model  covariance  matrix,  a  single 
point-source  signal  is  added.  This  is  done  by  sweeping  the  signal  bearing  parameters 
through  360°  in  azimuth  and  90°  in  elevation,  while  keeping  the  noise  power  fixed  and 
optimizing  the  signal  power  (Section  3.2).  Several  potential  signal  parameter  solutions 
corresponding  to  higher  values  of  the  cost  function  L  are  selected.  For  each  of  these 
candidate  solutions,  a  modified  gradient  ascent  technique  is  employed  to  jointly  refine 
both  the  noise  power  and  all  the  signal  parameters  including  spread  (see  Section  4.0). 
The  refined  parameter  set  yielding  the  highest  value  of  the  cost  function  is  then  retained 
and  all  other  refined  candidate  sets  are  thrown  out.  This  completes  stage  one. 

In  the  next  stage,  using  the  new  model  covariance  matrix  corresponding  to  the  refined 
model  parameter  set  from  the  previous  iteration,  a  new  set  of  point-source  signal  candi¬ 
dates  are  identified  using  the  azimuth  and  elevation  search  method  in  the  same  manner 
as  in  the  first  stage.  This  is  followed  by  another  joint  refinement  of  all  signal  and  noise 
parameters,  and  again,  only  the  refined  parameter  set  yielding  the  highest  value  of  L  is 
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retained.  This  completes  stage  two. 

The  process  described  for  stage  two  is  repeated  again  and  again  until  the  last  ( Mth ) 
signal  has  been  introduced  and  the  model  parameters  have  been  refined.  The  resultant 
model  parameters  are  then  taken  to  be  the  optimum  maximum  likelihood  estimates,  the 
determination  of  which  will  be  presented  later. 

An  extra  step  is  added  to  the  joint  refinement  procedure  when  the  bearing  of  a  candi¬ 
date  point-source  signal  falls  within  the  angular  spread  region  of  one  or  more  of  the  refined 
signals  selected  in  the  previous  stage.  In  this  case,  the  refinement  procedure  involving  the 
candidate  point-source  signal  is  repeated  with  the  following  special  modifications  to  the 
parameters  of  the  affected  signals  (including  the  candidate  signal  and  the  refined  signals 
with  which  it  “collides”,  but  no  others):  the  azimuth  and  elevation  spread  values  are  set 
to  zero,  and  the  signal  powers  set  equal  to  the  average  power  of  the  colliding  signals.  If 
this  modified  refinement  yields  the  highest  value  of  L,  then  the  new  refined  parameter  set 
is  retained,  otherwise  it  is  rejected. 

Performing  the  azimuth-elevation  search  at  each  step  using  a  point-source  signal  tends 
to  favour  the  detection  of  point-source  signals  over  spread-source  signals.  However,  since 
point-source  signals  are  of  the  most  interest  (i.e.,  for  high  latitude  HF  DF  applications,  a 
point-source  signal  is  more  likely  to  yield  the  true  bearing  of  the  transmitter  than  a  real 
spread-source  signal),  this  enhancement  is  not  a  bad  thing. 

The  identification  of  several  candidate  model  parameter  sets  during  the  search  and 
the  ensuing  refinement  to  identify  the  best  solution  is  a  modification  to  the  original  APM 
method  in  order  to  improve  the  probability  of  finding  the  correct  solution. 

The  modified  gradient  ascent  technique  used  for  parameter  refinement  is  also  a  de¬ 
viation  from  the  original  APM  method.  It  was  used  because  it  was  found  to  converge 
faster  than  both  the  technique  proposed  for  the  APM  method  and  the  standard  gradient 
ascent /descent  method. 

The  next  few  sections  provide  a  detailed  discussion  of  the  initial  noise  estimate,  the 
azimuth-elevation  search  procedure  for  the  introduction  of  new  signals,  and  the  modified 
gradient  ascent  technique  used  to  refine  the  signal  and  noise  parameters. 
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3.1  Initial  Noise  Estimate 


The  initial  noise  estimate  is  made  by  assuming  that  there  are  no  signals  present  and 
then  finding  the  value  of  a2  which  maximizes  the  cost  function  L.  The  no  signal  assump¬ 
tion  is  equivalent  to  assuming 

C  =  <r2CV  (32) 

Substituting  this  into  the  expression  for  cost  function  (10)  gives 

L0  =  —  ln(det(<72C,j))  -  trace(R(cr2C^)_1) 

=  -iVln(cr2)  -  In(detC^)  -  ^trace(RC“1)  (33) 

where  Lq  denotes  the  cost  function  value  for  zero  assumed  signals  and  N  is  the  number 
of  sensors. 

The  cost  function  can  be  maximized  (or  minimized)  by  differentiating  L0  with  respect 
to  cr2  and  then  finding  the  value  of  a2  which  sets  the  differentiated  expression  to  zero. 
Hence  differentiating  yields 

=  +  ^trace(RC^)  (34) 

where  u  —  a2,  and  then  setting  the  result  to  zero  and  solving  gives 

o2  =  -^:trace(RC71).  (35) 

For  this  value  of  a2  the  value  of  the  second  derivative  d2Lo / dv2  is 

-TV3 

trace (RC"1)2  <  °  ^ 

indicating  that  value  of  o2  given  by  (35)  produces  a  maximum  (not  minimum)  in  the  cost 
function  Lq. 

The  expression  for  the  optimum  value  of  a2  (in  the  no  signal  case)  can  be  inserted  back 
into  the  expression  for  the  cost  function  L0  (10)  to  give  a  new  expression  independent  of 
cr2  (but  where  the  optimum  value  of  a2  is  implicitly  assumed).  The  result  is 

Lq  =  N\n(N)  -  N  -  ln(det  C,)  -  ATln(trace(RC71)).  (37) 
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3.2  Azimuth-Elevation  Search 


In  the  azimuth-elevation  search,  a  new  point-source  signal  is  added  to  the  model 
covariance  matrix.  In  this  section,  the  search  procedure  is  developed  so  that  the  azimuth 
and  elevation  angles  of  the  new  signal  can  be  adjusted  without  having  to  explicitly  adjust 
the  signal  power  as  well.  Additionally,  although  the  model  noise  covariance  is  required 
in  this  search,  an  initial  estimate  of  the  noise  power  is  not. 


Beginning  with  the  known  model  covariance  represented  by  Ca,  the  addition  of  a  new 
point-source  signal  yields 

C  =  Ca  +  sm6mem  (38) 

where  the  indexing  begins  with  m  =  1  and  is  incremented  after  each  stage  until  m  =  M. 
The  inverse  of  the  matrix  C  can  be  found  using  the  augmented  matrix  inversion  lemma 
to  give 


0—1  _  ^-i—l  _  ^m6mCQ, 

°  l  +  ^egCr'e 


(39) 

where  it  is  also  assumed  here,  and  throughout  the  rest  of  this  report,  that  C  and  Ca 
are  positive  definite  so  that  C-1  and  C'1  exist  and  e^C~1em  >  0.  For  simplicity  a  new 
vector 

Qm  =  CQ  em 

is  introduced  so  that  the  expression  for  C-1  becomes 


(40) 


Q-l  _  0-1  _  smQmQ- 


H 


1  Sm&mQm 


(41) 


Plugging  the  two  expressions  for  C  and  C  1  into  (10)  and  simplifying 


H 

I'm  =  ln(det(CQ  +  smeme^))  -  trace(R(C”1  -  -  )) 


“1“  smGmqm 


=  -  ln(det(I  +  smqme^)  det(CQ))  -  trace(RCa1)  +  trace  ( \ 

\i  +  smemqm  J 

=  -  In((l  +  sme*qm)  det(Ca))  -  trace(RC;x)  +  )  (42) 

\l  +  sme"qro/ 

where  Lm  is  the  cost  function  after  the  mth  signal  is  introduced  into  the  model. 

To  reduce  the  search  to  only  the  azimuth  and  elevation  angles,  the  cost  function  can 
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be  differentiated  with  respect  to  the  signal  power  and  then  solved  to  find  the  optimum 
value  of  the  signal  power.  Differentiation  with  respect  to  sm  produces 


dLm 

dSm 


^mQrn  _  SmGm Q m Qm 

1  d"  f  d~  i  (1  d" 


and  setting  this  result  to  zero  then  solving  for  sm  gives 


Sm  — 


qmR<lm  -  emQm 


The  second  derivative  <92Lm/<9s^  for  this  value  of  sm  is 


-(e^qm)4 

(q£Rqm)2 


(43) 


(44) 


(45) 


indicating  that  this  value  of  sm  is  the  desired  maximum  value,  and  not  a  minimum. 


Finally,  the  expression  for  sm  can  be  reinserted  into  (42),  and  after  some  manipulation 
the  result  is  given  by 


-  ln(det  Ca)  -  trace  (RC~X)  -  ln(^Rqm)  +  q^Rq” 


em(lm 


e^q, 


-  1 


mHm 


j  ,q^Rqm,  q^Rq„ 


emQ m 


e£q m 


-  1 


(46) 


which  is  the  cost  function  evaluated  for  the  optimum  signal  power.  Note  that  in  the  above 
expression,  Lm-\  is  the  maximum  cost  function  value  from  the  previous  stage. 

Using  (46),  the  search  for  the  optimum  signal  parameters  can  be  carried  out  by  varying 
the  azimuth  and  elevation  angles  only.  The  maximization  process  can  be  simplified  even 
more  by  noticing  that  the  constant  terms  in  (46)  have  no  effect  on  the  maximization,  and 
that  maximizing  g(x)  -  lnp(rr)  with  respect  to  x  is  the  same  as  maximizing  g(x),  where 
g(x)  is  a  positive  real-valued  function.  Hence  maximizing  Lm  is  the  same  as  maximizing 


or 


^m) 


qmRqm 

emqm 


(47) 

(48) 
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The  function  Sm((j)m,  ipm)  is  not  only  the  specialized  cost  function  for  the  single  signal 
case  (or  additive  signal  case),  but  it  can  also  be  interpreted  as  the  whitened  beamformer 
spectrum.  To  see  this,  consider  that  in  estimation  problems  a  common  practice  is  to 
filter  the  data  in  such  a  way  as  to  whiten  the  noise  spectrum  (this  assumes  that  the 
noise  spectrum  is  known  a  priori).  Once  this  has  been  done,  optimal  estimators  designed 
for  white  noise  can  be  used.  This  simplifies  the  problem  of  having  to  design  different 
estimators  for  white  noise  and  coloured  noise. 

For  the  estimator  problem  at  hand,  CQ  can  be  treated  as  the  known  “noise”  covariance, 
and  assuming  it  to  be  full  rank,  then  there  will  exist  an  appropriate  N  x  N  spatial 
whitening  filter  W  such  that 

WWff  =  Ca-1  (49) 

(how  W  is  actually  determined  is  unimportant  here).  Using  this  filter,  the  whitened  data 
is  formed  according  to 

Xw  =  WHX  (50) 

which  leads  to  the  whitened  covariance  matrix 

Ru,  =  WffRW.  (51) 

To  account  for  the  affect  of  the  whitening  filter  on  the  response  of  the  sensor  array,  the 
steering  vector  e  is  modified  in  the  same  way  as  the  data  and  then  normalized  to  give 

WHe  _  WFe 

Gw  "  vPwPe  ~ 

The  whitened  beamformer  spectrum  is  then  given  by 

which  becomes  (48)  if  R^,  and  are  expanded  in  terms  of  (50)-(52). 

From  this  discussion,  it  can  be  seen  that  the  introduction  of  a  new  point-source  signal 
into  the  covariance  model  can  be  accomplished  using  a  simple  beamformer  to  search  in 
azimuth  and  elevation.  Signal  power  does  not  have  to  be  explicitly  tested.  Additionally, 
the  first  signal  can  also  be  determined  without  knowledge  of  the  noise  power  since  the 
maximum  of  Sm((f>m,  VVn)  remains  unchanged  whether  Ca  =  cr2Cv  is  used  or  Ca  =  Cv  is 
used. 


(52) 


(53) 


18 


4.0  PARAMETER  REFINEMENT 


In  the  previous  section,  initialization  of  the  model  parameters  was  discussed.  In  this 
section,  the  gradient  ascent  method  is  introduced  and  developed  to  further  refine  and 
optimize  the  model  parameter  estimates.  The  principle  strategy  employed  is  to  incremen¬ 
tally  adjust  the  values  of  the  model  parameters  according  to  the  local  gradient  of  the  cost 
function  until  a  maximum  in  the  cost  function  is  found  (i.e.,  climbs  the  “hill”  until  the 
peak  is  reached).  Since  this  approach  finds  the  local  maximum,  care  must  be  taken  in  the 
initial  model  parameter  estimates  to  ensure  that  the  local  maximum  corresponds  to  the 
global  maximum. 

Using  a  to  represent  any  of  the  real-valued  model  parameters  (e.g.,  a2,  or  sm ,  or  <f)m, 
etc.),  then  the  gradient  ascent  parameter  adjustment  at  each  step  i  is  performed  according 
to 

oii—i  4-  fj,aG(oii- 1)  (54) 

where  G(ai)  is  the  slope  or  gradient  with  respect  to  the  model  parameter,  and  jj,Q  is 
a  positive  real  value  used  to  control  the  step  size.  The  above  expression  is  performed 
simultaneously  for  all  the  adjustable  model  parameters  before  proceeding  to  the  next 
step,  that  is,  model  parameters  for  step  i  are  computed  using  only  the  model  parameters 
from  step  i  —  1. 

The  parameter  adjustment  step  size  ^  is  also  modified  at  each  step  (in  this  discussion 
the  subscript  a  is  used  only  when  discussing  the  step  size  for  a  single  model  parameter,  and 
dropped  when  discussing  the  step  sizes  for  all  the  model  parameters).  In  each  step,  after 
all  the  model  parameters  have  been  adjusted,  the  new  parameter  set  is  tested  to  see  if  the 
cost  function  is  improved.  If  L,  >  the  parameter  adjustment  is  considered  successful 
and  the  step  size  for  all  parameters  is  slightly  increased  (//  -4  1.2/z).  If  L;  <  Li+i,  the 
adjustment  is  considered  unsuccessful,  the  new  parameter  set  is  discarded,  and  the  step 
size  for  all  parameters  is  decreased  (/i  — »  fjt/3). 

In  addition  to  these  modifications  to  the  step  size,  a  special  modification  is  also  used  to 
improve  the  behaviour  of  the  gradient  ascent  technique  when  various  kinds  of  structures 
in  the  cost  function  are  encountered.  For  these  structures,  a  few  parameters  will  be  close 
to  their  optimum  values  but  their  corresponding  slopes  on  both  sides  of  these  optimum 
values  will  be  large,  while  for  other  parameters  their  optimum  values  are  further  away 
and  the  slopes  are  smaller  (e.g.,  consider  a  long  ridge  which  slowly  rises  along  the  length 
of  the  ridge  but  has  very  steep  sides).  The  result  is  that  the  closer  parameter  values  will 
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dominate  the  cost  function  (the  larger  slope  means  a  greater  change  in  the  cost  function 
for  a  given  increment),  and  the  values  of  /j,a  will  need  to  be  reduced  accordingly  to  keep 
these  parameter  values  in  the  vicinity  of  their  optimum  values.  With  a  small  value  of  /ia, 
the  parameters  that  are  further  away  take  a  large  number  of  steps  before  they  come  close 
to  their  optimum  values. 

One  way  to  detect  this  slow  convergence  condition  is  to  examine  the  behaviour  of  the 
gradient  function.  For  a  parameter  close  to  its  optimum  value,  the  gradient  will  oscillate 
between  being  positive  and  negative  as  the  adjustments  overshoot  back  and  forth  across 
the  optimum  value  (G{ai_l)G(ai_2)  <  0).  If  this  condition  is  detected  after  the  new 
parameter  set  has  been  deemed  successful,  then  the  step  size  jxa  is  replaced  by  fia/ 3. 


The  definition  of  the  gradient  function  depends  on  the  model  parameter  being  adjusted, 
however  before  getting  into  the  finer  details,  some  generalizations  can  be  made.  The 
gradient  of  the  cost  function  (10)  is  given  by 


G(a)  = 


dL 

da 


dln(detC)  dtrace(RC  *) 
da  da 


(55) 


This  can  be  simplified  somewhat  by  considering  each  of  the  two  terms  in  the  above 
expression  individually,  as  will  be  seen. 


Beginning  with  the  first 
covariance  matrix 


term,  consider  the  following  eigen-decomposition  of  the  model 

C  =  J2  Aiv*vf  (56) 

i=  1 


where  A,  represents  the  eigenvalues  and  v*  represents  the  eigenvectors  of  C  with  the 
properties  vf  v,  =  1  and  vfv*  =  0  for  i  ^  k.  Since  the  determinant  of  a  matrix  is  equal 
to  the  product  of  it’s  eigenvalues,  the  first  term  in  (55)  becomes 


d  ln(detC)  =  dE^M^i) 
da  da 


Performing  the  indicated  differentiation  and  utilizing  the  eigenvectors  and  their  properties 


jr 


20 


to  expand  and  rearrange  the  expression,  then 


51n(det  C) 
da 


N  d\- 

- 

=  f):A  . 

^  1  da  da 


&V^Vi 

where  — ^ — -  =  0  since  =  1 

5a  1 


“  1 1  + ^ 


=  traoe(E  £  ylKW»(^tv»  +  vk^y«  +  v,A»^f )) 


2=1  Jfc=l 


trace(£  v^vf  £  A»vf  +  vt^vf  +  v*A*M) 


=  trace(£vfA 

»=1  Jk=l 


trace(C  1-^— ). 


The  second  term  can  be  simplified  by  first  differentiating  the  equality  CC-1  =  I, 
which  gives 

dCCT1 

-&r=0  <59) 

where  0  is  a  matrix  of  zeros  having  the  same  dimensions  as  C.  Expanding  the  lefthand 
side 

dc  ,  ^ac-1 

+c^-=0  (6°) 

and  rearranging 

dC"1  ,dC  . 

"ST  =  'c  Sc  '  <61> 

Using  this  relationship,  and  carrying  out  the  indicated  differentiation  of  the  second  term 


in  (55),  then 


<9trace(RC  x) 
da 


dC~\ 

=  traceR— — ) 
da 

dC 

=  —  trace  (RC-1-r—  C-1) 
da 

=  — trace(C-1RC-1^^). 


Finally,  a  general  expression  for  the  gradient  can  be  found  by  substituting  (58)  and 
(62)  back  into  (57),  which  gives 

a/-* 

G(a)  =  -trace(C-1— )  +  trace(C_1RC_1— ) 


=  trace  ^(C_1R  —  I)C-1^-^  . 


Detailed  derivations  of  the  gradient  as  a  function  of  noise  power,  signal  power,  signal 
bearing,  and  bearing  spread  are  given  in  the  following  sections.  To  simplify  the  expressions 
that  result,  the  definition 

P  =  (C-1R  -  ^C"1  (64) 

is  used.  The  gradient  is  then  given  by 

dC 

G(a)  =  trace(P— ).  (65) 

da 

Additionally,  the  definition  for  the  model  covariance  C  used  in  the  following  derivations 
is  given  by  (31)  and  repeated  here  as 

M  Mg 

C  =  <7  C„  +  ^  Sm  y]  Pmi^Tnij  (66) 

m— 1  i—  1 

4.1  Noise  Power  Gradient 


Letting  a  =  a2  and  differentiating  the  model  covariance  matrix  expression  with  respect 
to  the  noise  power  gives 

=  C„,  (67) 
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and  incorporating  this  in  (65),  the  noise  power  gradient  becomes 

G(a2)  =  trace  (PC,,). 


4.2  Signal  Power  Gradient 

Differentiating  the  model  covariance  matrix  expression  with  respect  to  the  power  of 
the  mth  signal  gives, 

dC 

cT  (69) 

and  using  a  =  sm  in  (65),  then  the  signal  power  gradient  becomes 


G(sm)  =  trace(P^/3m(^mi,^mi)emie^). 

i—1 


4.3  Bearing  Gradient 

Differentiating  the  model  covariance  matrix  expression  with  respect  to  the  azimuth 
bearing  angle  of  the  mth  signal  gives 

d<t>m  m^m{  mu*mt){d<f>m  mi  +  emid(/>J-  (71) 

Using  the  definition  for  the  steering  vector  (15)  and  the  fact  that  differs  from  <f>m  by 
a  predetermined  constant  offset,  then 


where  Am,  is  a  diagonal  matrix  defined  by 


x0  cos  <j)mi  cos  ipmi  -  yo  sin  <f)mi  cos  Ipmi 
/ .  x  _  .  2tt  cos  (j)mi  cos  ipmi  -  yx  sin  <pmi  cos  ipmi 

—  j 

_  XN-1  COS  <f>mi  COS  1pmi  -  yN_x  sin  <j)mi  cos  1pmi 
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Using  these  definitions  and  oc  —  (j)m  in  (65),  the  azimuth  bearing  gradient  becomes 

Mg 

=  smtrace(P  ^  ’®,mi)(Amiemiemi  emjemiAm,-)).  (74) 

i= 1 

Similarly,  the  elevation  bearing  gradient  can  be  derived  to  get 

Mg 

=  smtrace(P  ^  mil  ^mi) (Umiem,emj  emiemiBmi))  (75) 

i= 1 

where  Bm,  is  a  diagonal  matrix  defined  by 

x0  sin  4>mi  sin  ipmi  +  y0  cos  <j)mi  sin  ^mi 

2t r  xi  sin  <f)mi  sin  ipmi  +  yi  cos  (j)mi  sin  tpmi 

diag(Bmi)  =  -j—  . 

£jv-i  sin  4>mi  sin  cos  $mi  sin 

A  special  modification  to  the  elevation  bearing  gradient  occurs  if  part  of  the  spread  region 
is  cut-off  by  the  horizon  or  90°  vertical.  The  expression  is  then  represented  by 

G(ipm)  =  smtrace  (p  ^  ^(^»^)em,eg 

y  1=1  Otym 

Mg  \ 

4”  P  )  .  ^ mi) (Pmi&mi&m-;  1  •  (77) 

Since  the  complete  derivation  requires  the  exact  definition  of  the  taper  function,  this 
derivation  is  discussed  separately  in  Section  4.5. 

4.4  Spread  Gradient 

Differentiating  the  model  covariance  with  respect  to  the  azimuth  or  elevation  spread 
of  the  mth  signal  yields 

dC  _  ^0/U*mi,*mi)  H 

dAm  Smh\  dAm 

+  +  Pm($mi,Vmi)emi^J  (78) 
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where  Am  represents  either  A^m  or  A^m.  Defining  Cmj  as  the  model  covariance  matrix 
for  each  grid  signal,  or 

Cmi  —  ^mi)®mi®mt5  (79) 


then  the  derivative  of  the  model  covariance  matrix  can  also  be  represented  in  a  more 
compact  form  as 


dC  ^dCmi 

dAm  Sm  dAm  ' 

The  corresponding  form  for  the  gradient  is  given  by 


(80) 


Ma 


'l9  . 

G{ Am)  =  smtrace(P  -^±)- 

i=l  0lXrn 


(81) 


Without  further  defining  the  grid  layout  and  taper  functions  it  is  not  possible  to  go 
any  farther,  so  for  this  report,  two  different  choices  of  grid  shape/taper  functions  were 
investigated.  They  were  as  follows: 


1.  A  rectangular  shaped  grid  using  uniform  grid  spacing  except  on  the  outside  edge, 
and  a  uniform  power  distribution  for  the  individual  point-source  signals. 

2.  An  elliptical  shaped  grid  with  uniform  grid  spacing  and  a  raised-cosine  power  dis¬ 
tribution  for  the  individual  point-source  signals. 

Although  the  second  choice  was  expected  to  more  closely  resemble  the  real  world,  assessing 
both  choices  provided  a  means  of  comparing  performance  when  very  different  patch  shapes 
were  assumed  for  the  same  data. 

In  the  following,  the  spread  gradient  is  defined  for  both  taper  functions  under  the 
assumption  that  neither  the  estimated  signal  elevation  bearing  nor  spread  angles  result  in 
grid  signals  below  the  horizon  or  above  vertical.  The  special  case  when  this  assumption 
is  not  true  is  considered  in  the  Section  4.5. 


4.4.1  The  Uniform  Rectangular  Grid 

Figure  2  shows  both  the  general  shape  of  a  rectangular  grid  and  the  arrangement  of 
point-source  signals  in  one  row  of  the  grid  (the  arrangement  in  a  column  is  the  same).  With 
the  exception  of  the  two  outside  signals,  the  grid  signal  locations  are  uniformly  spaced 
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with  uniform  power  and  the  central  signal  located  at  (4>m,  ^m)-  The  outside  signals  have 
their  power  levels  scaled  according  to  the  angular  distance  to  the  neighbouring  signal  in 
order  to  provide  a  smooth  transition  when  the  spread  is  varied  to  the  point  that  new 
rows/columns  of  point-source  signals  must  be  added  or  removed  from  the  grid. 

The  determination  of  the  values  of  the  taper  function  for  the  various  grid  signals  begins 
with  the  calculation  of  the  numbers  of  these  signals  as  a  function  of  the  spread  parameters. 
Along  each  row  (defined  by  constant  elevation  with  increasing  azimuth  angles  measured 
left  to  right),  but  excluding  the  outermost  or  nonuniformly  spaced  grid  signals,  there  will 
be 

M^=  2int(|fe-)  +  l  (82) 

grid  positions  where  d ^  is  the  angular  grid  spacing  within  the  row  (it  is  also  the  column 
spacing).  Similarly,  along  each  column  (defined  by  constant  azimuth  with  increasing 
elevation  measured  bottom  to  top),  there  will  be 

JK*.=2int(At)  +  l  (83) 

grid  positions,  where  c^m  is  the  spacing  within  the  column  (it  is  also  the  row  spacing). 
Ideally  d<pm  and  d7prn  are  chosen  to  be  as  small  as  possible,  but  since  this  increases  the  num¬ 
ber  of  computations  required,  reasonable  choices  are  one  fifth  the  array  3  dB  beamwidth 
(which  is  usually  different  in  the  azimuth  and  the  elevation  directions).  The  spacing 
requirement  is  discussed  and  evaluated  in  more  detail  in  Section  6.1. 

The  M^m  x  M^m  grid  signals  in  the  interior  of  the  grid  will  have  constant  power,  and 
the  corresponding  value  of  the  taper  function  can  be  defined  as  fom •  The  value  of  the 
taper  function  for  the  grid  signals  at  the  edges  is  a  function  of  the  spacing  at  the  edge 
relative  to  the  spacing  of  the  interior  signals.  Along  the  rows,  the  edge  spacing  is  given 
by 

da  —  —  —  1  )d(f>m).  (84) 

Similarly,  along  the  columns  the  edge  spacing  is  given  by 

de  =  —  —  1  )di/,m).  (85) 

Excluding  the  corner  signals  of  the  grid,  the  corresponding  value  of  the  taper  function  for 


26 


(b) 

Figure  2:  Geometry  of  the  uniform  rectangular  grid  showing  (a)  the  general  shape,  and 
(b)  an  example  of  the  position  and  power  of  the  point-source  signals  making  up  one  row 
of  a  grid. 
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the  edges  of  each  row  is  given  by 


Am  =  /Wp-  =  %(§*=■ -iU*m  +  1)  (86) 

Hm  *  Hm 

and  of  each  column  is  given  by 

Am  =  Am "7”"  =  -  M*.  +  1).  (87) 

For  the  four  corners  of  the  grid,  the  values  of  the  taper  function  will  be  a  function  of  both 
the  row  and  column  edge  spacing  values  leading  to 


p3m  —  fior 


dii 


%(^  -  Mtm  +  1)(^=-  -  Mim  +  1). 


Since  the  sum  of  the  values  of  the  taper  function  for  all  the  grid  positions  is  equal  to 
one,  then 

M (j)m  M1pm  Pom  +  +  2M(f>mP2m  +  4/?3  m  =  1.  (89) 

Substituting  for  /3lm,  /32m,  and  /?3m  using  (86)- (88),  and  then  rearranging  in  terms  of  /30m 
the  expression  becomes 

Am  =  + 1) 

V  d4>m 

+  Mtm(^  -  Mfm  +  1)  +  (4*=.  -  +  1)(^  -  Mfm  +  1)) 

d1pm  d<t>m  d”4,m  ) 

_  f  ^4>m  |  |  j  ^ 

\  d(pm  dxpm  c?^m  d^m 


Having  now  defined  the  shape  of  the  grid  and  the  taper  function,  it  is  now  possible  to 
continue  the  derivation  of  the  gradient  from  (78).  For  signals  in  the  interior  of  the  grid, 
only  the  power  varies  as  a  function  of  the  azimuth  or  elevation  spread,  so  that  dCmi/dAm 


becomes 


9Cmi  dP0m  jj 
dAm  ~  dAm  mi  mi 


(91) 


where  /?0m  is  the  value  of  (30m  calculated  for  the  mth  grid  signal. 


Using  the  result  from 
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(90)  and  carrying  out  the  derivative  operation  then 


dCr 


A) m 

d<t>m  d^m 


+  1  )emiemi- 


Excluding  the  corner  signals,  the  signals  on  the  top  and  bottom  edges  of  the  grid  also 
vary  only  in  power  as  the  azimuth  spread  varies.  Noting  that  the  power  is  given  by  p2m, 
which  can  be  expanded  in  terms  of  p0m  using  (87),  then  for  the  top  and  bottom  edge 
signals  the  appropriate  expression  is 


dCmt  _  ,  1\.  aH 

r)A  A  'A  '  ^v®****®!?** 

C'A0m  “0m  “0m 

which  is  similar  to  (92). 


For  signals  on  the  left  and  right  edges  of  the  grid  (and  again  excluding  the  corner 
signals),  both  the  power  and  signal  direction  vary  according  to  the  azimuth  spread.  The 
relationship  between  azimuth  spread  and  signal  direction  for  the  right  and  left  edge  signals 
is  given  by 

<f>mi  =  <t>m  ±  (94) 

where  “+”  is  used  for  the  right  edge  and  is  used  for  the  left  edge.  Using  this  relation¬ 
ship  plus  the  fact  that  the  power  level  is  given  by  /?lm,  and  partially  differentiating  Cmi 
with  respect  to  A^m,  then 


dCr 


dAd 


OP}™.  e  .e#.  ±  0lm  ( A  .g  -H 
dA^  2  ^’n"m‘l'Srni'£rni 


■eH  A  ■) 


where  A  was  previously  defined  in  Section  4.3.  Using  the  definition  for  /?lm  given  in  (86) 
and  the  definition  for  p0m  given  in  (90)  in  order  to  carry  out  the  rest  of  the  differentiation, 
the  above  result  becomes 


dC 


mi 


dAd 


Por, 
2  cL 


1  -  2 'Puni-?*-  +  1) 

.  “0m  / 


e  eH- 


0U 


(A, 


■eFA 


i).  (96) 


Finally,  for  the  signals  at  the  corners  of  the  grid,  the  differentiation  proceeds  in  a 
similar  manner  except  that  the  power  level  is  given  by  p3m  which  is  defined  in  terms  of 
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A>m  by  (88).  The  resultant  equation  is  given  by 


dCmi 

d^m 


—  /Am  A)  A  i  A*2L(  A  e  eH  -  e  eHA  ■)  (97') 

—  I  n  ,  ,  V  j  ^  A/  I  emiemt  ^  o  •  J 

\2d4>m  Hm  di>m  J  2 


where,  as  before,  “+”  is  used  for  the  corners  on  the  right  edge  of  the  grid  and  is  used 
for  corners  on  the  left  edge  of  the  grid. 

Using  the  same  procedure  to  differentiate  the  covariance  matrix  with  respect  to  the 
elevation  spread,  then  for  signals  located  in  the  interior  of  the  grid 

_  ~AL  /  &4>m  ■  -|\  h  /qo\ 

Excluding  the  grid  signals  at  the  corners,  then  for  signals  on  the  right  and  left  edges 


dCmi 


Am  Am  /  <Pm  .  1  \  H 

j  \  j  >  L)^mt^rnit 

Hm  a4>m 


and  for  signals  on  the  top  and  bottom  edges 

1?=-  =  |r-  f1  -  2/M^  +  1))  ±  %(Bm(em(e"  -  em(e"  Bmj)  (100) 

2d, V  J  2 


where 


V’mi  =  1pm  ± 


(101) 


and  where  “+”  is  used  for  the  top  edge  and  is  used  for  the  bottom  edge.  Finally,  for 
the  grid  signals  at  the  corners 


dCmi 


@0m.@3m  | 

dipm 


emieTni 


±  ^(B 


mi&mi&rni  ^ 


mi^mi 


H  R  \ 

(102) 


Having  differentiated  the  covariance  matrix  in  terms  of  all  the  component  signals  of 
the  grid,  the  gradient  G(A<pm)  is  computed  using  (81)  where  dCmi/dA^m  is  defined  by 
one  of  equations  (92),  (93),  (96),  or  (97)  depending  on  the  location  of  the  signal  in  the 
grid. 


Similarly,  the  gradient  G(A^m)  is  also  computed  using  (81)  where  dCmi/dA^m  is 
computed  using  one  of  equations  (98),  (99),  (100),  or  (102)  depending  on  the  location  of 
the  signal  in  the  grid. 
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4.4.2  The  Raised-Cosine  Elliptical  Grid 


Figure  3  shows  both  the  general  shape  of  an  elliptical  grid  using  a  raised-cosine  taper 
function  and  the  arrangement  of  point-source  signals  in  one  row  of  the  grid  (the  arrange¬ 
ment  in  a  column  is  the  same).  With  the  taper  function  going  to  zero  at  the  edges  of  the 
region,  no  special  edge  spacings  are  required  to  introduce  or  remove  new  signals.  However, 
for  smaller  spread  values  it  becomes  necessary  to  reduce  the  grid  spacings  so  that  at  least 
three  point-source  signals  are  used  in  the  center  row  or  center  column  of  the  grid.  Other¬ 
wise  the  resultant  signal  could  end  up  being  composed  of  a  signal  point-source  signal  in 
each  row  and/or  column  with  no  way  to  compute  the  spread  gradient.  Additionally,  for 
signals  where  the  spread  region  drops  below  the  horizon,  or  above  vertical,  modifications 
are  required  which  are  discussed  in  Section  4.5. 

For  convenience,  it  will  be  initially  assumed  that  both  the  azimuth  and  elevation  spread 
values  (A^  and  A^)  are  large  enough  that  the  spacing  parameters  and  can  be  chosen 
to  be  some  constant  fraction  of  the  horizontal  and  vertical  beamwidths  of  the  antenna 
array.  The  case  involving  smaller  spread  values  will  then  be  dealt  with  afterwards. 


The  taper  function  for  the  elliptical  grid  is  defined  by 


Pmi  — 


Cm(l  +  COs(27!Tmi))  if  Tmi  <  0.5 


otherwise 


(103) 


where  /3mi  is  used  as' a  shorthand  way  of  writing  /3m($mi,  ^mi),  cm  is  a  normalization 
factor  so  that  (3mi  =  1  and  can  be  determined  from 


E*Si(!  +  cos(27rrmfc)) 


(104) 


and  rmi  is  the  normalized  distance  from  the  center  of  the  grid  expressed  mathematically 


rrni  =  V^mi  +  ^mi 


0m  0mi  \  .  ( 0m  0mi  \ 


(105) 


To  avoid  unnecessary  computations,  the  grid  is  only  generated  in  the  elliptical  region 
defined  by  rmi  <  The  resultant  elliptical  grid  then  measures  A^m  along  the  azimuth 
axis  from  left  to  right  and  A^m  along  the  elevation  axis  from  top  to  bottom  (or  slightly 
less  than  these  dimensions  depending  on  the  spacing  parameters).  Additionally,  if  the 
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(a) 


Q. 


^  4>+2d* 

Azimuth  Angle 

(b) 

Figure  3:  Geometry  of  the  raised-cosine  elliptical  grid  showing  (a)  the  general  shape, 
and  (b)  an  example  of  the  position  and  power  of  the  point-source  signals  making  up  one 
row  of  a  grid. 
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grid  extends  below  0°  or  above  90°  in  elevation  (which  is  an  unrealistic  condition  for  the 
kinds  of  signals  being  dealt  with  in  this  report),  these  parts  of  the  grid  are  removed. 


Based  on  the  fixed  spacing  for  all  the  signals  used  in  the  grid,  the  intermediate  result 
(91)  found  for  the  rectangular  grid  is  also  appropriate  for  the  elliptical  grid  and,  after 
making  the  appropriate  modifications,  is  shown  here  as 


dCr 


dp, 


"mi  h 


dAm  8Am 

Using  the  definition  for  pmi  for  the  elliptical  grid,  and  letting  Am  =  A^m  then 
dpmi  47 r2cm 


(106) 


Ma 


dAT 


A4>n 


$mi  sinc(2rmj)  -  pmi  £  $2mk  sinc(2 rmk) 


k-\ 


(107) 


and  the  azimuth  spread  gradient  (106)  becomes 


dCr, 

dA * 


47T2Cm 

A4>m 


Ma 


S> 


mi 


sinc(2rmi)  Ani  y]  sinc(2rmfc) 

k= 1 


e  e11 


(108) 


where  sinc(x)  =  sin(jra)/(jra)  and  sinc(0)  =  1  .  Letting  Am  =  Afci  the  equivalent 
elevation  spread  gradient  becomes 


dCmi  47 T2Cm  (  M9  \ 

=  ~A~  sinc^2r^')  “  P™  £  *lk  sinc(2rmfc)  emie"  .  (109) 

Wvn  \  fc-1  J 

These  last  two  expressions  can  then  be  substituted  back  into  (81)  to  compute  the  gradients 
G(A4>m)  and  G(A^m),  respectively. 

If  the  azimuth  and/or  elevation  spread  value  is  reduced  sufficiently,  the  above  approach 
fails  because  the  number  of  rows  and/or  columns  in  the  grid  is  reduced  to  one.  For 
example,  if  the  spread  A $  <  2d^  in  Figure  3,  then  there  would  be  a  single  point-source 
signal  at  (j)  in  each  row  of  the  grid.  At  this  point,  then,  the  corresponding  value  of  pmi  will 
not  change  as  the  given  spread  is  further  reduced  so  that  the  corresponding  differential 
becomes  zero.  This  condition  occurs  when 


and/or 


A4>m 


(110) 

(111) 
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where  c^m  and  d^m  are  fixed  values. 

To  overcome  this  difficulty,  modeling  of  the  spread  signal  is  changed  so  that  the  central 
row  and  column  contain  at  least  three  point-source  signals.  This  is  done  by  using  the 
adjustable  azimuth  grid  spacing 

4.  =  for  <  44.  <112) 

instead  of  the  fixed  spacing  c^m,  and  using  the  adjustable  elevation  grid  spacing 

^Ipm  ^V’ro  ^  ^dipm  (H3) 

instead  of  the  fixed  spacing  d .  An  example  for  the  case  of  narrow  spreading  in  azimuth 
is  shown  in  Figure  4.  The  conditions  placed  on  e^m  and  d'^m  in  the  preceding  expressions 
have  been  chosen  to  ensure  continuity  between  the  narrow  spread  model  and  the  wide 
spread  model  since  for  =  4e^m  and/or  A^m  =  4d^m,  the  two  models  are  identical. 


Azimuth  Angle 

Figure  4:  Position  of  point-source  signals  in  a  raised-cosine  elliptical  grid  used  for  narrow 
spread  signals. 

For  the  narrow  spread  model,  the  angles  of  the  outer  point-source  signals  are  always  in 
the  same  relative  position  with  respect  to  the  central  point-source  and  the  spread  angle, 
so  that  the  corresponding  values  of  (3mi  become  constant.  For  example,  assuming  that 
A^m  <  4 d,),  then  the  values  of  (f>mi  are  chosen  from  the  set 

tm,  K  +  (114) 

Since  4>mi  =  {<j)m  -  ^mi)/A^m,  then  the  corresponding  values  of  are  chosen  from  the 
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set 


(115) 


,1  1 .. 

{  4’  °’  4 

which  is  independent  of  A^m .  Based  on  the  relationship  between  <&mi  and  (5mi  expressed  in 
(103)-(105),  then  /?mi-  will  also  be  independent  of  the  azimuth  spread.  Similarly,  assuming 
that  A^m  <  4d<p  then  the  values  of  (f>mi  are  chosen  from  the  set 


{*»  -  o,  i>m  - 


(116) 


Since  =  (ipm  —  ipmi)/A^m,  the  resultant  values  of  are  also  given  by  (115),  and  as 
before,  this  leads  to  f3mi  being  independent  of  the  elevation  spread. 

Deriving  the  gradient  for  the  narrow  spread  case  is  very  similar  to  deriving  the  bearing 
gradient  since  only  the  point-source  signal  directions  are  affected  by  changes  in  spread, 
and  not  the  value  of  the  taper  function  f3mi.  Hence  using  the  results  from  the  bearing 
gradient  case,  then  for  A^m  <  Ad^  the  azimuth  spread  gradient  is  then  given  by 


dC 

dAs„ 


M„ 


—  ® m  ^  '  fimi  ( 


de 


■mi  H 


em,-  +  e„ 


de11- 

ucmi 


2=1 


a  a.  ^  ™  a  Ad 


) 


(117) 


where 


(  4  Amiem,  for  —  4 


de. 


dAd 


=  0 


for  <S>mi  =  0 
for  $77ii  =  t 


For  Aipm  <  Ad$  the  elevation  spread  gradient  is  given  by 


(118) 
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(119) 


where 


f  4  for  ^mi  —  4 


vCmi 


dA 


=  < 


Ipn 


lB  e 


for  $mi  =  0 


for  $mj;  =  7 


(120) 


v  4  mi  mz  ~  4 

The  matrices  A mi-  and  Bmi  were  defined  previously  by  (73)  and  (76),  respectively,  in 
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Section  4.3. 


4.5  Out-of-Bounds  Signals 

This  section  deals  with  the  problem  of  modeling  the  spread  region  with  a  planar  array, 
when  part  of  the  grid  is  adjusted  so  that  it  falls  below  the  horizon  or  above  vertical  (out- 
of-bounds).  Without  compensation,  the  out-of-bound  signals  are  folded  back  in-bounds 
corrupting  the  shape  of  the  signal  model  spread  region  (i.e.,  grid  signals  below  the  horizon 
act  as  if  the  elevation  angle  is  while  signals  above  the  vertical  act  as  if  the  elevation 

angle  is  V’mi  —  90°.  The  compensation  required  depends  on  the  taper  function  used. 


4.5.1  The  Uniform  Rectangular  Grid 

For  the  uniform  rectangular  taper  function,  the  problem  of  out-of-bound  grid  signals 
can  be  easily  handled  by  readjusting  the  elevation  angle  and  spread  parameters.  For 
spreading  below  the  horizon  (tpm  —  A^m/2  <  0°)  the  modifications  are 

Vi  =  iWw  +  %)  (121) 

Ai,  =  2<P'm  (122) 

where  the  ip'm  and  A^m  are  the  updated  values,  and  for  spreading  above  vertical  (ipm  — 
A^m/2  >  90°)  the  modifications  are 

Vi  =  |(*»  -  %■)  +  45°  (123) 

=  2(90° -C)-  (124) 

These  readjustments  avoid  the  need  to  derive  expression  (77)  since  the  readjustment  can 
be  done  without  affecting  the  shape  or  distribution  of  the  spread  region. 
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4.5.2  The  Raised-Cosine  Elliptical  Grid 


For  the  raised-cosine  elliptical  taper  function  things  are  not  as  easily  handled  as  for 
the  uniform  rectangular  taper  function  due  to  its  nonuniform  nature.  Instead  the  part  of 
the  grid  that  is  out-of-bounds  is  first  removed.  However,  this  leads  to  discontinuities  in 
the  gradient  functions  when  grid  signals  are  added  or  deleted  in  this  manner.  To  handle 
this  transition,  a  new  row  of  grid  signals  are  introduced  right  at  the  boundary  as  shown 
in  Figure  5.  The  power  of  these  boundary  signals  are  determined  in  much  the  same  way 
as  done  for  the  edge  signals  of  the  uniform  rectangular  grid.  The  model  covariance  matrix 
is  then  formed  using 


M  / 

sm  I  flmjemjGmj  *b 
m=l  \  j 

where  the  shorthand  form  f3mi  is  again  used  for  $my),  and  the  summation  with 

respect  to  the  index  j  is  assumed  to  include  all  boundary  grid  signals  while  the  summation 
with  respect  to  the  index  k  is  assumed  to  include  all  other  grid  signals  which  are  within 
the  boundaries.  To  remain  consistent  throughout  the  rest  of  this  section,  the  indices  j  and 
j'  will  be  used  to  denote  a  parameter  associated  with  the  grid  signals  on  the  boundary, 
and  k  and  k'  will  be  used  to  denote  a  parameter  associated  with  grid  signals  inside  the 
boundaries.  It  is  also  assumed  that  the  elevation  spread  region  is  narrow  enough  that  the 
removed  grid  rows  are  either  all  below  the  horizon  or  all  above  vertical,  but  not  both. 


C  —  a2  Cv  +  ^2 


Pmk&mkemk 


(125) 


Elevation  Angle 


Figure  5:  Position  of  point-source  signals  in  a  raised-cosine  elliptical  grid  where  the 
spread  region  goes  below  the  horizon.  Dotted  arrows  show  angular  positions  of  rows  that 
have  been  removed. 


The  value  of  the  taper  function  for  An*  is  computed  in  the  same  way  as  before,  namely, 

An*  =  ^(1  +  cos  2tt  rmk)  (126) 
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while  the  adjusted  value  of  the  taper  function  f)mj  is  computed  using 


a. 


Pmj  =  Cm(  1  +  -j“)(l  +  COS  2*  pmj) 


where  dm  is  the  grid  elevation  spacing  parameter  defined  as 

(  4m  for  A^m  >  44, 


dm  —  \ 


for  A <  Ad , 


Aprr 


(127) 


(128) 


and  am  is  the  angular  elevation  difference  between  the  boundary  and  the  first  deleted  grid 
row  as  shown  in  Figure  5.  More  specifically,  it  is  defined  according  to 
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VVl 


90°  -  ipn 


below  the  horizon 


above  vertical 


(129) 


where  ipmv  is  defined  as  the  elevation  angle  of  the  removed  (or  virtual)  row  of  grid  signals 
closest  too,  but  below  the  horizon  or  above  vertical  (the  corresponding  elevation  angles 
il>mj  are  either  0°  or  90°).  The  normalizing  constant  then  becomes 


Cm  — 


53(1  +  COS  27ITmfc)  +  (1 +  ^)  53(1+ COS  2?r Pmj,) 

k  «m  j> 


(130) 


The  definition  for  rm*  is  given  by  (105)  and  the  definition  of  pmj  is  the  same  as  rm*  except 
that  the  elevation  angle  used  in  this  case  is  ipmv  (instead  of  ipmj),  so  the  definition  becomes 


pmj  —  — 


(131) 


Additionally,  T mj  is  used  in  place  of  in  the  above  expression  to  indicate  the  modifi¬ 
cation  to  the  elevation  angle. 


Using  these  new  definitions,  the  calculation  of  the  noise,  signal  power,  azimuth  bearing, 
and  azimuth  spread  gradients  proceed  in  the  same  manner  as  before.  The  gradients  for 
the  elevation  bearing  and  elevation  spread,  however,  need  modification. 

Starting  with  the  elevation  bearing,  and  given  the  modifications  made  to  the  signal 
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grid,  the  gradient  given  by  (77)  can  be  expressed  as 


G(ipm)  =  smtrace[P^^^emie^. 

j  Otym 


+  p£ 
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OPmk 
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^mk^mk  "b  Pmk{^mk^mk^mk  emk&mk^mk) 


(132) 


Performing  the  derivative  operations  indicated  then 
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(134) 


where  given  a  choice,  “+”  is  used  for  below  horizon  conditions  and  ”  is  used  for  above 
vertical  conditions.  The  elevation  gradient  expressed  in  (132)  can  now  be  computed  by 
substituting  (133)  and  (134)  for  the  derivatives. 

Before  continuing  on  to  the  spread  gradients,  it  is  useful  to  distinguish  between  the 
two  separate  cases,  namely  modeling  signals  with  a  wide  elevation  spread  where 


A*»>4d*m,  (135) 

and  modeling  signals  with  a  narrow  elevation  spread  where 

<  4^.  (136) 
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The  main  difference  is  that  in  the  wide  case  the  grid  elevation  spacing  is  fixed  while  in 
the  narrow  case  the  grid  elevation  spacing  changes  as  a  function  of  the  elevation  spread 
angle  (see  (128)).  Since  these  differences  lead  to  different  derivations,  the  two  cases  are 
handled  separately. 

Starting  with  the  wide  elevation  spread  case,  the  gradient  is  given  by 


—  smtrace  P )  1 


VmjVmj 


dA 


+  pE^|.  (137) 


/ 


The  derivatives  of  (3mj  and  pmk  are  practically  identical  to  the  corresponding  derivative 
for  (3mi  in  the  previous  section  (107)  except  for  minor  modifications  to  account  for  the 
boundary  signals.  With  these  modifications,  then 
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(139) 


The  gradient  for  the  azimuth  spread  (wide  case)  can  then  be  computed  by  substituting 
(138)  and  (139)  for  the  appropriate  derivatives  in  (137). 

For  the  narrow  azimuth  spread  case  the  gradient  is  given  by 


G{A^m)  =  smtrace[P^3 


dfimj  H 


dA,emjemj 


j 


+p£ 


dflrnk  H  ,  o  /  ^ernk  H  .  „  ®emk 
^rnk^mk  T  Pmk  l  ^  emfc  *r  emk^X 


k 


dA,, 


(140) 


where  demk/dA^m  is  defined  by  (120).  Using  the  fact  that  ipmj  =  ipm  -  dL  for  below 
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horizon  signals,  ipmj  =  ipm  +  d!^m  for  above  vertical  signals,  and  d!^m  =  A^m  /4  then 


A nj  =  Cm  4^-(l  +  COS  2lTrmj) 


and 


Cm  — 


4/4, 


Ed  +  cos  2nrmk)  +  ir^Ed  +  cos  2nrmji ) 

k 


(141) 


(142) 


where 


f  A 


Vm=  { 


below  the  horizon 


(143) 


[  90°  —  'ipm  above  vertical 


The  differentiation  of  (3mj  then  proceeds  as 


d(3mj 

dA^Pm 


dcm  .  4  vm  .  N 

A  (1  “I"  cos  27irmj)  cm  g  (1  +  cos  27rrmj) 


16k 


A 


— J^(l  +  cos  2nrmj)  £(1  +  cos  2?rrm;v)  -  -2JL 

7V 


ftmj 


E& 


A 


mj 


A. 


In  a  similar  manner,  the  differentiation  of  Anfc  yields 


dPmk 

d^ipm 


Pmk 


j' 


(144) 


(145) 


Based  on  the  derivations  here  and  those  in  the  previous  section,  the  gradient  for  the 
azimuth  spread  (narrow  case)  can  be  computed  from  (140)  by  substituting  (120),  (144), 
and  (145)  as  appropriate. 
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5.0  SML  ALGORITHM  SUMMARY 


The  following  two  sections  summarize  the  main  SML  routine  and  the  gradient  ascent 
routine.  The  gradient  ascent  routine  was  summarized  separately  owing  to  its  complexity, 
otherwise  the  resultant  combined  routine  would  have  been  more  difficult  to  follow. 

In  the  main  SML  routine,  two  different  constant  values  are  introduced  which  affect 
algorithm  performance.  The  first  is  a  30%  threshold  level  used  to  determine  whether  a 
maximum  value  found  for  the  function  tpm)  is  significant  with  respect  to  the  global 

maximum  value.  A  higher  value  could  be  used  to  reduce  the  number  of  maximum  values 
which  exceed  this  threshold,  however  through  numerous  simulations  it  has  been  found 
that  this  reduces  accuracy  for  conditions  where  angular  separations  between  signals  are  a 
beamwidth  or  less. 

The  second  constant  value  is  the  grid  spacing  which  is  chosen  as  30%  of  the  array 
beamwidth.  The  reasons  for  this  choice  are  discussed  later  in  Section  6.1. 

5.1  Main  SML  Routine 

1.  Compute  the  data  covariance  matrix  using  (11) 

R  =  4 XXH .  (146) 

K 

2.  Set  the  stage  counter  index:  m  =  0. 

3.  Estimate  the  normalized  noise  covariance  matrix  based  on  using  theoretical  con¬ 
siderations  or  practical  measurements.  For  example,  the  theoretical  covariance  for 
white  Gaussian  noise  is  found  according  to  (18) 

c„  =  (147) 

Alternatively,  if  measurements  of  noise  only  data  are  available  then  the  normalized 
noise  covariance  matrix  is  given  by  (19) 

XXg 

"  trac  e(XXH)' 


42 


4.  For  the  moment,  assume  the  data  contains  no  signals  and  estimate  the  noise  power 
using  (35) 

°2  =  ^trace(R<V)>  (149) 

then  generate  the  initial  model  covariance  matrix  using  (32)  with  C  =  Co  giving 

C0  =  ct2C„.  (150) 


5.  Increment  the  index:  m  =  m  +  1. 

6.  Create  a  list  of  candidate  azimuth  and  elevation  bearings  (within  the  angular  ranges 
of  interest)  yielding  local  maximum  values  within  30%  of  the  global  maximum  of 
the  function  defined  by  (48)  where  Ca  =  Cm_i,  or 


<Sm(^mj  ^m) 


-’-i 


(151) 


The  matrix  Cm_i  represents  the  model  covariance  matrix  generated  at  stage  m  —  1, 
and  the  array  response  or  steering  vector  is  given  by  (15) 


g 3  X  (X°  sin  0m  COS  1pm  +yo  COS  <pm  COS  1pm  ) 
gj'xfc1  sin0m  COS  Ipm+Vl  COS  <pm  COS  1pm) 


(152) 


7.  For  each  list  entry  containing  candidate  signal  bearings,  include  the  corresponding 
signal  power  calculated  using(40),  (44),  and  (47),  or  equivalently 


1^'m(0m)'0m)  1 

pffO-1  p 


(153) 


and  the  spread  values 


^ <Pm  ^Ipm  0. 


(154) 


Also  include  all  the  signal  and  noise  model  parameters  from  the  optimum  model 
computed  at  stage  m  —  1.  Hence  each  list  entry  will  contain  the  initial  estimates  for 
a  possible  solution  involving  one  new  signal,  m  -  1  previously  established  signals, 
and  noise. 


8.  If  for  any  list  entry  the  new  candidate 


signal  bearing  coincides  (collides)  with  the 
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scattering  region  of  any  of  the  established  signals,  perform  the  following: 

(a)  create  a  new  list  entry  by  repeating  all  the  parameters  from  the  colliding  entry; 

(b)  reset  the  signal  power  of  all  colliding  signals  in  the  new  entry  to  the  average 
signal  power  of  these  same  signals;  and 

(c)  set  the  azimuth  and  elevation  spread  values  of  all  colliding  signals  to  zero. 

9.  Refine  the  noise  and  signal  parameter  values  of  each  list  entry  using  the  gradient 
ascent  technique.  If  the  uniform  rectangular  taper  function  is  employed,  make 
adjustments  to  the  appropriate  parameters  if  out-of-bound  conditions  are  detected 
as  outlined  in  Section  4.5.1. 

10.  Update  the  model  covariance  matrix  Cm  associated  with  each  list  entry  according 
to  the  current  model  parameter  refinements  using  either  (31) 

m 

Cm  =  o2Cn  +  'Esn'Epn($ni,*1u)enie&  (155) 

n= 1  i= 1 

where  the  values  of  /?n(4>ni,  \kni)  are  computed  as  discussed  in  Section  4.4,  or  if  the 
raised-cosine  elliptical  taper  function  is  employed  and  out-of-bounds  conditions  are 
detected  then  use  (125) 


m  ( 

—  G  @t)  -f-  ^  )  Sn  I  ^  '  PnjGnjGnj  +  )  ^  Pnk&nk ^nk 

n=l  \  j  k 

where  /?„_,■  and  /?„*  are  computed  as  discussed  in  Section  4.5.2. 

11.  Also  compute  the  cost  function  value  of  the  refined  model  for  each  list  entry  using 
(10)  with  C  =  Cm  which  yields 


Lm  =  -  ln(det  Cm)  -  trace  (RCmx).  (157) 


Note  that  if  the  number  of  list  entries  is  excessive,  pruning  of  this  list  to  improve 
computational  speed  can  be  accomplished  by  performing  Steps  9-11  using  only  20 
iterations  of  the  gradient  ascent  technique  and  then  retaining  the  list  entries  cor¬ 
responding  to  the  four  highest  cost  values.  Steps  9-11  are  then  repeated  for  the 
four  remaining  entries  using  the  full  number  of  iterations  for  the  gradient  ascent 
technique. 
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12.  Select  the  model  parameters  in  the  list  corresponding  to  the  highest  cost  function 
value.  These  are  considered  the  optimum  model  parameters  estimates  for  stage  m. 

13.  Repeat  Steps  5-13  as  long  as  m  <  M. 

14.  Output  the  model  parameter  results  from  stage  M. 


5.2  Gradient  Ascent  Routine 


1.  Save  the  value  of  the  cost  function  for  the  current  unrefined  model:  Lmax  =  Lm. 


2.  Initialize  the  gradient  ascent  loop  parameters  according  to  the  following 


loop  =  0 
/A?  =  £ 

Mst  =  £ 

/%  =  £ 
Aty*  =  £ 
=  f 
=  ? 


where  k  =  1, m  and  £  is  the  initial  step  size.  A  suggested  approach  to  determine 
an  appropriate  initial  step  size  is  to  compute  the  initial  gradients  of  the  azimuth  and 
elevation  bearings  and  the  initial  gradients  of  the  bearing  spreads  for  all  signals  in  the 
current  model  using  Steps  4-5.  The  maximum  absolute  value  gmax  is  selected  from 
these  gradients  and  the  initial  step  size  is  then  computed  using  £  =  0.05<f)bw/gmax 
where  <j)bw  is  the  3  dB  azimuth  beamwidth  of  the  antenna  array.  This  restricts  the 
initial  adjustment  of  the  bearing  and  spread  parameters  to  no  more  than  5%  of  the 
azimuth  beamwidth. 


3.  Increment  the  loop  counter:  loop  =  loop  +  1. 

4.  Compute  the  quantity  (64) 


P  =  (C-1R-I)C-1. 


(158) 
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5.  Compute  the  gradients 


G(a 2)  =  trace  (PC^)  (159) 


Mg 

G(sk)  =  trace(P53^fc($fei,^i)efcie^)  (160) 

1=1 

Mg 

G{<t>k)  =  SktraceiPY.M^ki^kdi^ki^kieki-ehie^Aki))  (161) 

1=1 

Mg 

G(ipk)  =  sktra.ce(P'%2/3k($kh'&ki)(BkiekieEi- ekiefjBki))  (162) 

1=1 

(from  (68),  (70),  (74),  and  (75),  respectively)  for  k  =  1  and  where  A ki  and 
Bfcj  are  diagonal  matrices  defined  by  (73)  and  (76)  respectively,  which  are  repeated 
here  as 


diag(Afci) 


x0  cos  4>ki  cos  ifiki  ~  Vo  sin  (j>ki  cos  ipki 
X\  cos  <f>ki  cos  ipki  -  yi  sin  <j)ki  cos  ipki 


Xpf— i  COS  (pki  COS  Ipki  VN—  1  Sin  <pki  COS  'ipki 


(163) 


diag(Bjfci) 


x0  sin  (pki  sin  ipki  +  yo  cos  4>ki  sin  ipki 
X\  sin  (pki  sin  ipki  +  y\  cos  <pki  sin  ipki 


xN-i  sin  <pki  sin  ipki  +  2/v-i  cos  (pki  sin  ipki 


(164) 


Since  the  gradients  for  G( and  G(A1pk)  cannot  be  easily  summarized  in  one 
line,  refer  to  Section  4.4.  Additionally,  for  out-of-bounds  conditions  (which  affects 
G(ipk)  and  G(A^k)),  refer  to  Section  4.5. 


6.  Update  the  model  parameters  by  adapting  (54)  as  appropriate  and  letting  i  =  loop 
to  get 


<72 

->  cr2  +  pvG(cr2) 

(165) 

sk 

— >  sk  +  pSkG(sk) 

(166) 

(pk 

-»  (pk  +  p^>kG((pk) 

(167) 
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fa  — >  fa  +  ^kG(fa) 

Afa  -»  A,pk  +  fi^kG(A^k) 

Aipk  ->  A^  +//a^G(A^J 


(168) 

(169) 

(170) 


for  fc  =  1, ...,  ra.  If  the  uniform-rectangular  taper  function  is  used  and  out-of-bounds 
conditions  are  detected,  adjust  the  appropriate  model  parameters  as  outlined  in 
Section  4.5.1. 


7.  Compute  the  model  covariance  matrix  using  (31) 

m 

Cm  =  a2Cv  +  52  Sn  5Z  Pn($ni,  ^ni)^ni^ni  (171) 

n=l  i~ 1 

where  the  values  of  A»($n»j^nt)  and  the  corresponding  grid  of  point-source  sig¬ 
nals  are  computed  as  discussed  in  Section  4.4  (see  equations  (82)-(90)  for  a  rect¬ 
angular  grid  and  equations  (103)-(105)  for  an  elliptical  grid)  Alternatively,  if  the 
raised-cosine  elliptical  taper  function  is  employed  and  out-of-bounds  conditions  are 
detected  then  using  (125) 


m  / 

Cm  =  cr2Crj  ~H  Sn  I 
n=l  \ 

where  finj,  /?„*  and  the  corresponding  grid  of  point-source  signals  are  computed 
as  discussed  in  Section  4.5.2  (see  equations  (126)-(131)).  The  grid  spacing  values 
are  chosen  as  fam  =  O.Zfaw  and  =  0.3^,  where  faw  is  the  3  dB  elevation 
beamwidth  of  the  antenna  array. 

8.  Update  the  cost  function  using  (10) 


H 

njc";cnj 


“b  Pnk^nk^nk 
k 


(172) 


Lm  =  -  ln(det  Cm)  -  trace(RCm1).  (173) 


9.  If  the  new  cost  value  is  less  than  the  old  cost  value  ( Lm  <  Lmax)  then: 

(a)  reset  the  model  parameters  and  the  model  covariance  matrix  to  their  original 
values  at  the  beginning  of  the  loop; 

(b)  reduce  all  the  step  size  parameters  by  a  factor  of  3;  and 

(c)  if  loop  <  200  go  to  Step  3,  otherwise  continue  on  to  Step  10. 
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Otherwise  if  ( Lm  >  Lmax)  then: 


(a)  save  the  new  cost  value:  Lmax  =  Lm; 

(b)  increase  all  step  size  values  by  a  factor  of  1.2; 

(c)  compare  all  gradient  values  to  the  corresponding  value  from  the  previous  loop 
and  if  there  has  been  a  sign  change  (negative  to  positive  or  vice  versa)  decrease 
the  associated  step  size  value  by  a  factor  of  3;  and 

(d)  if  any  of  the  azimuth  or  elevation  bearings  changed  more  than  5  x  10-6  degrees 
in  Step  8  and  loop  <  200,  then  go  to  Step  3,  otherwise  continue  on  to  Step  10. 

10.  If  loop  <  5  ,  then  reinitialize  the  step  size  parameters  (using  the  same  method  as 
in  Step  2)  and  go  to  step  3.  Otherwise  continue  on  to  the  next  step. 

11.  Output  the  final  model  parameter  results. 
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6.0  SML  CONTROL  PARAMETER  SETTINGS 


In  the  previous  sections,  the  SML  algorithm  has  been  described  in  some  detail.  A  few 
of  the  more  practical  considerations  have  not  yet  been  discussed  adequately.  This  includes 
the  control  parameter  settings  for  the  grid  spacing  to  be  used  for  the  signal  model,  the 
amount  of  data  that  should  be  collected  to  form  the  data  covariance  matrix,  the  choice 
of  the  taper  function,  and  the  assumed  number  of  signals. 

To  investigate  these  control  parameters,  numerical  simulations  were  carried  out  to 
determine  their  most  appropriate  settings  or  range  of  settings.  In  the  simulations,  the 
signals  received  at  the  antenna  array  consisted  of  a  single  point-source  signal,  representing 
a  signal  from  the  desired  great  circle  direction,  plus  two  stronger  spread  signals,  represent¬ 
ing  signals  off  the  great  circle  bearing  scattered  by  closely  following  ionospheric  patches. 
The  spread  signals  were  modeled  using  a  large  number  of  point-source  signals  spread  with 
random  amplitudes,  phases,  and  bearings.  Unless  otherwise  specified,  the  spread  signals 
were  set  up  to  have  an  elliptical  shape  with  a  raised-cosine  power  distribution.  The  spread 
signals  were  also  uncorrelated  from  sample  to  sample  by  computing  every  sample  using 
different  random  values.  Ground  reflections  and  local  site  multipath  effects  were  not  in¬ 
cluded.  Noise  was  uncorrelated  (both  spatially  and  temporally)  white  Gaussian  noise. 
The  relevant  parameters,  as  measured  at  the  receiver,  are  given  in  Table  1. 


Table  1:  Simulation  Signal  Parameters 


Signal 

$ 

A* 

Power 

1 

50° 

10° 

0° 

0° 

-20  dB 

2 

O 

O 

00 

H 

30° 

30° 

10° 

0  dB 

3 

210° 

25° 

15° 

5° 

0  dB 

noise 

- 

- 

- 

- 

-15  dB 

The  antenna  array  geometry  was  as  shown  in  Figure  6.  This  geometry  was  investigated 
in  [12]  and  [13]  and  found  to  have  very  good  characteristics  for  direction  finding.  Assuming 
an  ideal  free  space  response,  for  a  given  signal  bearing  the  3  dB  azimuth  beamwidth  of 
the  array  is  relatively  constant  with  respect  to  the  azimuth  bearing  but  varies  with  the 
elevation  bearing  according  to 

<(>bw  =  7.8°/|cos^|  for  \ip\  <  90°  (174) 
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N  ✓ 


Figure  6:  Three  dimensional  view  of  the  antenna  array.  Each  grid  square  has  a  dimension 
of  1A  x  1A. 


Signal  Elevation  Angle  (degrees) 


Figure  7:  Antenna  array  elevation  beamwidth  as  a  function  of  elevation  comparing  the 
simulated  response  (solid  line)  with  the  sin-1  xp  predicted  response  (dashed  line). 
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for  the  array  size  shown.  Note  that  measuring  the  azimuth  beamwidth  at  or  near  ip  =  90° 
is  effectively  meaningless.  The  3  dB  elevation  beamwidth  is  given  by 

ipBW  =  7.8°/|sin^|  for  |-0|  >  30°.  (175) 

For  elevation  bearings  below  30°  the  3dB  beamwidth  is  somewhat  more  complicated  as 
shown  in  Figure  7.  The  failure  of  the  above  expression  at  the  lower  elevation  angles  is 
due  to  the  distortion  of  the  main  lobe  in  the  antenna  pattern  at  the  elevation  angle  ip 
by  the  reflection  of  this  lobe  at  the  elevation  angle  —ip.  For  example  a  2-dimensional 
x-y  array,  using  the  free  space  assumptions,  has  a  symmetrical  gain  pattern  for  elevation 
angles  above  and  below  the  horizon,  i.e.,  a  main  lobe  and  a  reflection  lobe.  At  low  signal 
elevations,  the  main  and  reflection  lobes  begin  to  merge.  Using  the  3  dB  criteria,  the 
lobes  are  considered  merged  when  the  minimum  gain  between  the  two  lobes  is  greater 
than  -3  dB  the  maximum  gain.  For  the  array  configuration  shown  this  occurs  at  21.2°. 
At  this  point,  the  beamwidth  effectively  doubles.  For  even  lower  signal  elevations  the 
merged  lobes  move  closer  together  so  the  beamwidth  actually  decreases. 

When  processing  the  simulated  data  using  the  SML  algorithm,  the  default  values  for 
the  parameters  under  investigation  (unless  otherwise  specified)  were:  a  grid  spacing  of 
0.3  times  the  3  dB  beamwidth  of  the  antenna  array  ( d $  =  0.3 (pBw  and  =  O.SipBw)', 
a  data  block  size  of  K  =  100  samples;  the  raised-cosine  elliptical  taper  function;  and  an 
assumed  number  of  signals  equal  to  the  actual  number  of  signals. 

The  processed  results  were  quantified  in  two  ways.  The  first  was  the  measurement  of 
the  failure  rate  of  signal  bearing  estimates,  and  the  second  was  the  measurement  of  the 
accuracy  of  the  successful  estimates.  In  this  report,  a  bearing  estimate  was  considered 
to  be  a  failure  if  it  deviated  from  the  true  signal  bearing  by  more  than  half  the  array 
beamwidth  (taking  into  account  both  the  azimuth  and  elevation  beamwidths).  Accu¬ 
racy  was  determined  by  calculating  the  root-mean-squared  (RMS)  error  of  the  estimates 
according  to 

RMS  Error  =  jj  YjjL  -  ^n)2  +  ftL  -  V’m)2  (176) 

where  the  summation  was  performed  for  all  H  successful  estimates  of  signal  m.  Other 
specifics  of  the  simulations  and  processing  were  dependent  on  the  aspect  of  the  SML 
algorithm  being  investigated,  and  are  discussed  in  the  following  three  sections. 
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6.1  Signal  Model  Grid  Spacing 

The  angular  spacing  between  point-source  signals  used  in  the  signal  model  grid  must 
be  sufficiently  narrow  for  the  model  to  adequately  represent  a  real  spread  signal.  However, 
making  this  spacing  too  narrow  can  unnecessarily  increase  the  number  of  computations 
involved  in  generating  these  models  (the  increase  in  computations  is  inversely  proportional 
to  the  square  of  the  spacing).  Hence  it  is  useful  to  determine  the  widest  spacing  that  can 
be  used  before  introducing  too  much  error  due  to  the  poorer  modeling  that  results. 

Since  the  resolving  power  of  any  antenna  array  is  a  function  of  the  beamwidth  of  the 
array,  a  natural  choice  would  be  to  use  a  spacing  value  which  is  some  fraction  of  the 
beamwidth.  To  evaluate  the  effect  of  changing  the  azimuth  and  elevation  grid  spacing 
values,  a  series  of  100  simulations  were  carried  out  generating  100  sets  of  data.  The  data 
set  was  then  processed  using  grid  spacings  which  were  varied  from  0.1  beamwidths  up  to 
1.0  beamwidths  in  0.1  beamwidth  increments  in  both  azimuth  and  elevation  relative  to 
the  appropriate  azimuth  and  elevation  beamwidth.  The  results  are  shown  in  Figure  8. 

In  generating  Figure  8  the  true  bearings  were  taken  to  be  those  estimated  using  the 
0.1  beamwidth  spacing  instead  of  the  actual  bearings  given  in  Table  1.  This  is  the  only 
occasion  that  it  was  done  this  way.  The  intent  was  to  find  the  widest  spacing  that  yielded 
the  same  results  (or  reasonably  close  to  the  same  results)  as  if  the  signal  model  had 
been  generated  using  infinitesimally  narrow  spacings.  A  spacing  of  0.1  beamwidths  was 
considered  sufficiently  narrow  for  these  purposes.  Note  also  that  perturbations  observed 
in  the  plotted  results  (relative  to  the  general  trends)  are  almost  certainly  statistical  effects 
introduced  by  using  finite  (versus  infinite)  data  sets. 

Examining  the  relative  failure  rate,  it  remains  low  up  until  0.8  beamwidths.  (Using 
the  actual  signal  bearings  from  Table  1,  instead  of  the  values  estimated  for  a  spacing 
of  0.1  beamwidths,  the  failure  rate  varies  from  18%  to  20%  for  the  same  beamwidths). 
Above  0.8  beamwidths  the  increase  is  more  dramatic.  From  these  results  it  would  appear 
that  spacings  of  up  to  0.8  beamwidths  could  be  used  with  reasonable  performance  still 
achieved.  For  the  rest  of  the  results  provided  in  this  report  a  more  conservative  spacing 
of  0.3  beamwidths  was  used. 
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(a) 


(b) 

Figure  8:  Effect  of  grid  spacing  on  SML  estimation  performance  relative  to  the  perfor¬ 
mance  for  a  grid  spacing  of  0.1  beamwidths.  The  plots  show  (a)  the  relative  failure  rate 
for  each  of  the  three  bearing  estimates,  and  (b)  the  corresponding  relative  accuracy  of  the 
successful  estimates.  The  solid  lines,  dashed  lines,  and  dash-dot  lines  correspond  to  the 
results  for  the  point-source  signal  at  <j>  —  50°,  the  spread-source  signal  at  <f>  =  180°,  and 
the  spread-source  signal  at  </)  =  210°,  respectively. 
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6.2  Data  Block  Size 


In  the  assessment  of  the  effects  of  ionospheric  scattering  patches,  it  has  been  found 
that  signals  arriving  from  one  part  of  a  scattering  region  appear  to  be  uncorrelated  with 
signals  arriving  from  any  other  part  of  the  scattering  region.  To  take  advantage  of  this 
requires  collecting  enough  samples  of  sensor  data  over  a  long  enough  time  period  that  the 
signals  become  reasonably  decorrelated.  Assuming  that  the  samples  are  taken  at  intervals 
which  exceed  the  correlation  time  constant  of  the  data,  then  the  question  becomes:  how 
many  samples  are  enough? 

To  answer  this  question,  a  series  of  simulation  data  sets  were  produced.  Each  set 
consisted  of  100  trials  with  a  given  number  of  sensor  samples  (or  block  size)  K  collected 
for  each  trial,  which  was  varied  from  set  to  set.  To  account  for  the  most  obvious  effects  of 
averaging  inherent  in  processing  larger  block  sizes  (which  leads  to  an  equivalent  increase 
or  processing  gain  of  10  log  AT  dB  in  the  SNR),  the  noise  level  was  also  adjusted  according 
to 

Noise  Level  —  10 log#  —  35  dB.  (177) 

This  noise  level  compensation  was  done  to  avoid  masking  other  effects  which  might  not 
otherwise  be  apparent.  The  results  are  shown  in  Figure  9. 

Inspection  of  the  failure  rate  reveals  that  it  is  high  for  very  small  block  sizes  (e.g., 
K  =  1, 2, 3)  but  decreases  rapidly  to  zero  as  the  block  size  increases.  The  accuracy  results 
are  interesting.  Ignoring  the  effects  of  simple  averaging,  then  for  the  point-source  signal, 
there  is  no  significant  change  in  accuracy  as  the  block  size  increases  since  it  can  be  well 
modeled,  even  for  a  single  sample.  Additionally  it  is  already  spatially  decorrelated  from 
the  other  two  signals  given  the  angular  separation  of  several  beamwidths.  The  bearing 
estimation  accuracy  for  the  spread  signals  decreases  at  a  relatively  constant  rate  with 
a  reduction  in  the  error  by  a  factor  of  10  as  the  block  size  is  increased  from  1  to  100. 
The  improvement  is  a  reflection  of  the  fact  that  the  larger  block  sizes  enhance  the  effects 
of  decorrelation  on  which  the  signal  modeling  relies.  Given  that  the  spread  signals  are 
ten  times  greater  in  power  than  the  point-source  signal,  this  improvement  would  only  be 
expected  to  continue  until  the  RMS  errors  were  on  the  order  of  ten  times  less  than  that 
of  the  point-source. 

This  upper  limit  on  accuracy  as  a  function  of  block  size  is  confirmed  in  the  simulation 
results  shown  in  Figure  10.  In  this  case  only  the  two  spread  signals  were  simulated  and 
the  noise  level  was  15  dB  higher  than  specified  in  (177).  Comparing  the  RMS  error  results 
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Figure  9:  Effect  of  the  number  of  data  samples  (block  size)  on  SML  estimation  perfor¬ 
mance  showing  (a)  the  failure  rate  for  each  of  the  three  bearing  estimates,  and  (b)  the 
relative  accuracy  of  the  successful  estimates.  The  solid  lines,  dashed  lines,  and  dash-dot 
lines  represent  the  point-source  signal  at  (f>  =  50°,  the  spread-source  signal  at  </>  =  180°, 
and  the  spread-source  signal  at  <f>  =  210°,  respectively.  Note  that  the  effect  of  simple 
averaging  as  a  function  of  block  size  has  been  removed. 
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Block  Size  (K) 

(b) 

Figure  10:  Effect  of  the  number  of  data  samples  (block  size)  on  the  estimation  of  the 
spread  signals  for  a  lower  SNR  (15  dB  lower  than  in  the  previous  example)  showing  (a) 
the  failure  rate  for  the  two  spread  signal  bearing  estimates,  and  (b)  the  relative  accuracy 
of  the  successful  estimates.  The  dashed  and  dash-dot  lines  represent  the  spread-source 
signal  at  <j)  =  180°  and  the  spread-source  signal  at  <j>  =  210°,  respectively.  Note  that  the 
effect  of  simple  averaging  as  a  function  of  block  size  has  been  removed. 
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to  those  in  Figure  9(a),  below  K  =  12  the  results  are  very  similar  despite  the  difference 
in  SNR.  Above  this,  the  RMS  error  plots  in  Figure  10  flattens  out  as  suggested  already, 
indicating  that  sufficient  decorrelation  has  occurred  for  the  spread  signals  to  be  accurately 
modeled.  For  these  higher  values  of  K,  the  dominant  error  mechanism  is  noise. 

The  total  effect  of  increasing  the  block  size  on  RMS  bearing  error  can  be  determined 
by  including  the  effects  of  simple  averaging.  The  effect  of  averaging  on  accuracy  is  found 
by  considering  the  equivalent  increase  in  SNR  and  then  utilizing  the  results  from  Section 
7.1  to  relate  this  to  the  RMS  bearing  error.  Combining  this  with  the  results  observed 
here,  then  the  block  size  can  be  related  to  the  RMS  error  e  according  to 


r  K 

for  small  K 

.  Vk 

for  large  K 

where  a  constant  SNR  is  assumed.  In  the  above  expression,  small  K  means  that  not 
enough  samples  are  available  for  sufficient  signal  decorrelation,  while  large  K  means  there 
are  enough  samples  and  noise  effects  dominate.  For  a  point-source  signal  the  delineation 
between  small  and  large  is  K  =  1. 

Based  on  the  failure  rate  results,  K  =  N  (number  of  samples  equal  to  the  number  of 
sensors)  is  a  reasonable  lower  limit  for  the  number  of  samples.  A  more  suitable  choice 
would  be  to  select  a  value  of  K  large  enough  that  sufficient  signal  decorrelation  has 
occurred,  although  this  would  involve  a  careful  assessment  of  the  expected  signal  envi¬ 
ronment  in  terms  of  the  expected  signal  SNR’s  and  spreading  values.  For  the  rest  of  the 
simulations  that  follow  in  this  report  K  =  100  was  deemed  suitable. 

One  final  observation  is  that  wider  spatial  spreading  of  the  received  signal  results  in 
poorer  RMS  accuracy  of  the  bearing  estimate  (everything  else  being  equal).  This  effect 
is  discussed  in  more  detail  in  Section  7.2  . 


6.3  Taper  Function 

The  choice  of  taper  function  is  ideally  based  on  the  shape  of  the  ionospheric  scattering 
region  and  the  density  of  signal  power  across  this  region.  However,  since  the  ionosphere  is 
highly  dynamic,  it  is  likely  that  there  will  be  differences  in  these  characteristics  between 
different  scattering  regions,  or  even  differences  for  the  same  region  observed  at  different 
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periods  of  time. 


To  investigate  the  effect  of  mismatches  between  the  assumed  signal  model  and  the 
actual  signal  spreading  distribution,  three  simulation  data  sets  were  produced,  each  con¬ 
sisting  of  100  trials.  The  first  data  set  was  produced  using  the  signal  parameters  listed  in 
Table  1  and  using  the  raised-cosine  elliptical  power  distribution  to  generate  the  simulated 
signals.  The  second  data  set  was  produced  using  the  signal  parameters  listed  below  in 
Table  2  and  used  the  uniform  rectangular  power  distribution  to  generate  the  simulated 
signals.  The  third  data  set  was  produced  also  using  the  signal  parameters  listed  below 
in  Table  2,  but  using  the  asymmetric  power  distribution  shown  in  Figure  11  to  generate 
the  simulated  signals.  For  all  simulations  the  SNR  was  set  to  -12  dB  (lowered  from  the 
default  value  of  -15  dB  to  increase  the  probability  of  failures). 


Table  2:  Simulation  Signal  Parameters  for  Uniform  Rectangular  and  Asymmetric  Spread 
Distributions 


Signal 

<t> 

ip 

A  <t> 

Power 

1 

50° 

10° 

0° 

0° 

-20  dB 

2 

180° 

30° 

18.9° 

6.3° 

0  dB 

3 

210° 

25° 

9.5° 

3.2° 

0  dB 

noise 

- 

- 

- 

- 

-12  dB 

To  simplify  the  comparison  of  results,  the  spread  values  in  Table  2  were  chosen  so  that 
processing  any  of  the  data  sets  using  the  same  SML  algorithm  variant  produces  nearly  the 
same  mean  spread  values.  For  example,  using  the  SML  algorithm  with  the  raised-cosine 
elliptical  taper  function  (denoted  SML(e))  produced  mean  spread  values  that  were  the 
same  as  those  listed  in  Table  1  for  all  three  data  sets.  This  means  that  the  spread  results 
were  correct  for  the  first  data  set  but  were  overestimated  for  the  other  two  data  sets 
the  results  due  to  the  modeling  mismatch.  Similarly,  using  the  SML  algorithm  with  the 
uniform  rectangular  taper  function  (denoted  SML(u))  produced  mean  spread  values  that 
were  the  same  as  those  listed  in  Table  2  for  all  three  data  sets.  In  this  case  the  spread 
results  for  the  first  data  set  were  underestimated  while  the  results  for  the  other  two  data 
sets  the  results  were  correct. 

To  process  the  three  data  sets,  both  the  SML(e)  and  SML(u)  algorithm  variant 
were  employed.  For  both  variants  grid  spacings  of  d<pm  =  0.1  <f>Bw  and  d^m  =  0.1  ipBw 
beamwidths  were  used  because  it  was  found  that  for  spacings  of  0.3  beamwidths  the  failure 
rate  significantly  increased  when  using  SML(u)  (but  not  SML(e)).  Table  3  shows  the  over¬ 
all  RMS  bearing  error  and  failure  results  for  each  signal.  The  results  are  shown  in  order, 
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Figure  11:  Asymmetric  power  distribution  used  for  spread  signal  simulations. 


that  is,  the  corresponding  bearing  angles  are  $  =  (50°,  180°,  210°)  and  ^  =  (10°,  30°,  25°). 


Table  3:  Effect  of  Taper  Function  on  SML  Performance 


Signal  Type 

SML(e)  Results 

SML(u)  Results 

RMS  Errors 

Failure  Rate 

RMS  Errors 

Failure  Rate 

Elliptical 

Uniform 

Asymmetric 

(0.50°,  0.20°,  0.14°) 
(0.41°,  0.17°,  0.13°) 
(0.42°,  0.95°,  0.53°) 

(1%,  0%,  0%) 
(2%,  0%,  0%) 
(1%,  0%,  0%) 

(0.54°,  0.20°,  0.15°) 
(0.42°,  0.16°,  0.12°) 
(0.44°,  0.77°,  0.47°) 

(5%,  0%,  0%) 
(2%,  0%,  0%) 
(2%,  0%,  0%) 

Examining  the  results  for  the  point-source  (at  0  =  50°),  SML(e)  performed  better 
than  SML(u)  regardless  of  the  spread  distribution  types  of  the  other  two  signals.  For 
the  spread  signals  (at  <f>  —  180°  and  210°) ,  the  SML  algorithm  employing  the  matched 
taper  function  performed  better  for  both  spread  signal  types  (i.e.,  SML(e)  performed  best 
for  signals  using  the  raised-cosine  elliptical  spread  distribution,  and  SML(u)  performed 
best  for  signals  using  the  uniform  rectangular  spread  distribution).  The  degradation 
that  resulted  from  using  an  improperly  matched  taper  function  was  almost  insignificant 
in  terms  of  accuracy  for  the  uniform  and  elliptical  signal  types  although  there  was  an 
increase  in  the  failure  rate  for  SML(u).  For  the  asymmetric  signal  type,  the  accuracy  in 
estimating  the  spread  signals  was  significantly  degraded  although  the  errors  were  still  less 
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than  1°  .  The  poorer  accuracy  was  due  to  the  signal  models  for  both  SML(e)  and  SML(u) 
being  biased  towards  the  peak  in  the  asymmetric  spatial  power  distribution  of  the  spread 
signals.  This  biasing  was  greater  for  SML(e)  than  SML(u).  Removing  the  biasing  yields 
error  results  very  similar  to  the  results  for  the  elliptical  and  uniform  signal  types. 

Based  on  preceding  results,  it  seems  clear  that  the  best  performance,  not  surprisingly, 
is  obtained  when  the  taper  function  is  matched  to  the  actual  signal  spread  distributions, 
although  the  SML  approach  seems  to  be  reasonably  robust  when  mismatches  occur  (at 
least  for  the  type  of  array  tested  in  this  report  -  a  larger  array  with  a  narrower  beamwidth 
could  lead  to  a  greater  sensitivity  to  modeling  errors).  Given  the  superior  performance  of 
the  raised-cosine  elliptical  taper  function  for  point-source  signals,  and  the  greater  likeli¬ 
hood  that  this  taper  function  would  be  better  matched  to  real  world  signals,  it  was  used 
exclusively  throughout  the  rest  of  this  report. 


6.4  Assumed  Number  of  Signals 

A  common  difficulty  with  many  advanced  direction  finding  algorithms  is  determining 
the  correct  number  of  signals  intercepted  in  the  receiver  pass  band.  In  this  report,  unless 
otherwise  stated,  it  is  assumed  that  the  number  of  signals  is  known.  Given  that  the  number 
of  signals  may  sometimes  be  incorrectly  chosen,  it  is  useful  to  know  what  effect  this  has 
on  performance  of  the  SML  algorithm.  Using  the  default  settings  for  the  simulations  (as 
discussed  in  Section  6.0)  to  produce  100  trials,  the  data  was  then  processed  assuming 
from  one  signal  up  to  six  signals.  The  overall  results  are  shown  in  Figure  12. 

Inspecting  these  results  reveals  several  things.  When  the  number  of  signals  is  under¬ 
estimated,  then  either  the  strongest  signals  are  detected  only,  or  two  or  more  signals  are 
wrongly  identified  as  being  a  single  signal.  For  example,  in  the  case  where  only  one  signal 
is  assumed,  the  weakest  signal  (the  point-source  signal  at  <j>  =  50°)  is  never  seen.  The 
spread  signal  at  (f>  —  210°  is  sometimes  correctly  identified  while  the  rest  of  the  time  the 
two  spread  signals  are  mistakenly  identified  as  a  single  signal  with  an  intermediate  bear¬ 
ing.  In  terms  of  initial  detection,  the  spread  signal  <p  =  210°  is  favoured  over  the  spread 
signal  at  <j>  =  180°  since,  even  though  the  two  signals  have  equal  power,  the  signal  at 
(f>  =  210°  has  less  spreading  giving  it  a  higher  power  density  (more  power  per  solid  angle). 
When  two  signals  are  assumed,  the  two  stronger  spread  signals  are  properly  estimated, 
but  the  weaker  point-source  signal  still  remains  undetected. 

When  the  number  of  signals  is  overestimated,  either  spurious  results  are  generated,  or 
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Figure  12:  Effect  of  the  assumed  number  of  signals  on  SML  estimation  performance 
showing  (a)  the  failure  rate  for  each  of  the  three  bearing  estimates,  and  (b)  the  relative 
accuracy  of  the  successful  estimates.  The  solid  lines,  dashed  lines,  and  dash-dot  lines 
represent  the  point-source  signal  at  <f>  =  50°,  the  spread-source  signal  at  <j>  =  180°,  and 
the  spread-source  signal  at  4>  =  210°,  respectively.  The  results  for  spurious  signals  are  not 
shown. 
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a  single  spread-source  signal  is  identified  as  two  spread-source  signals.  The  splitting  of 
the  spread-source  signals  explains  the  worsening  error  as  the  number  of  assumed  signals 
is  increased  past  three  (the  bearings  of  the  split  signals  tend  to  be  offset  on  either  side  of 
the  true  signal). 

In  general,  from  the  point  of  view  of  estimator  performance,  it  is  better  to  overestimate 
the  number  of  signals  than  underestimate  the  number  of  signals,  since  overestimating 
yields  all  the  signal  directions.  The  main  problem  with  overestimating  is  the  generation 
of  spurious  signal  directions  and  errors  due  to  splitting  of  the  spread  region. 
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7.0  SML  PERFORMANCE 


The  previous  section  dealt  with  modifying  various  control  parameters  and  the  effect 
these  modifications  had  on  the  performance  of  the  SML  algorithm.  In  this  section,  the 
effect  of  various  uncontrolled  parameters  and  their  effect  on  the  performance  of  the  SML 
algorithm  is  studied.  These  parameters  are  discussed  in  the  next  few  sections  and  include 
SNR,  signal  spreading,  and  angular  spacing  between  signals.  Where  appropriate,  the 
results  using  the  MUSIC  DF  algorithm  are  also  shown  for  comparison  purposes. 


7.1  Effect  of  Noise 

As  was  done  when  investigating  the  effects  of  various  control  parameters,  the  effect  of 
noise  was  investigated  through  simulation.  In  the  first  series  of  simulations,  a  single  point 
source  signal  at  (<f>,  ip)  =  (180°,  30°)  was  generated  and  the  SNR  was  varied  from  -20  dB 
to  +40  dB  in  2  dB  increments.  One  hundred  trials  were  generated  for  each  SNR  setting. 
In  the  second  series  of  simulations,  the  point-source  signal  was  changed  to  a  spread-source 
signal  with  spread  parameters  (A*,  A^)  =  (30°,  15°),  but  all  other  parameters  remained 
the  same.  The  results  are  shown  in  Figure  13. 

Several  features  of  the  results  are  worth  pointing  out.  The  SNR  at  which  the  failure 
rate  dramatically  departs  from  0%  is  called  the  threshold  SNR.  For  the  point-source  this 
occurs  between  -12  and  -10  dB  and  for  the  spread-source,  it  is  slightly  worse,  occurring 
between  -8  and  -10  dB.  The  poorer  performance  for  the  spread  signal  (including  both 
threshold  and  accuracy)  is  a  function  of  the  amount  of  spreading  and  is  discussed  in  more 
detail  in  the  next  section. 

The  scales  used  to  display  the  RMS  bearing  errors  in  Figure  13(b)  were  chosen  because 
they  linearize  the  accuracy  results  for  the  point-source  signal  above  the  threshold  SNR. 
In  this  region,  for  every  20  dB  increase  in  SNR  the  RMS  bearing  error  is  reduced  by  a 
factor  of  10.  Written  mathematically,  the  relationship  is  expressed  as 

e  oc  V  SNR  (179) 

where  e  is  the  RMS  bearing  error.  The  spread-source  also  exhibits  the  same  behaviour 
between  -8  and  0  dB,  but  above  this  SNR  the  error  begins  to  level  out  as  the  uncertainty 
in  the  signal  model  begins  to  dominate  the  error  (as  discussed  previously  in  Section  6.2). 
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(a) 


(b) 


Figure  13:  Effect  of  the  signal-to-noise  level  on  SML  estimation  performance  showing  (a) 
the  failure  rate,  and  (b)  the  accuracy  of  the  successful  estimates.  The  solid  lines  represent 
the  point-source  signal  and  the  dashed  lines  represent  the  spread-source  signal. 
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At  higher  SNR’s  then,  the  only  way  to  improve  accuracy  would  be  to  use  a  larger  block 
size  assuming  the  signal  was  well  modeled. 

To  determine  whether  including  the  azimuth  and  elevation  spread  parameters  in  the 
SML  signal  model  degrades  performance  when  estimating  point-source  signals,  the  point- 
source  simulation  results  were  reprocessed  using  the  MUSIC  algorithm.  The  results  were 
virtually  identical  to  the  SML  results  indicating  that  any  degradation  that  does  occur  is 
insignificant. 


7.2  Effect  of  Signal  Spreading 

The  effect  of  signal  spreading  was  investigated  through  simulation  by  measuring  ac¬ 
curacy  as  the  azimuth  and  elevation  spreads  of  the  generated  signal  was  varied.  In  the 
simulation  a  single  signal  at  ((f),  ip)  =  (180°,  30°)  was  used  with  three  different  SNR’s, 
namely,  0,  15,  and  30  dB.  The  azimuth  spread  was  varied  from  1°  to  30°  according  to 
the  sequence  =  1°,  2°,  4°,  6°,  8°,  10°,  15°,  20°,  25°,  30°  for  each  value  of  the  SNR.  The 
elevation  spread  was  set  to  one  half  the  azimuth  spread  (A^  =  A^/2).  One  hundred  trials 
were  generated  for  each  spread  and  SNR  setting.  The  results  are  shown  in  Figure  14. 


Figure  14:  Effect  of  signal  spread  on  SML  accuracy  for  SNR’s  of  0,  15,  and  30  dB. 
Elevation  spread  was  varied  according  to  one  half  the  azimuth  spread. 

The  effect  of  signal  spreading  on  accuracy  can  be  attributed  to  two  factors:  signal 
model  uncertainty  and  the  filter  effect.  The  first  factor,  signal  model  uncertainty  arises 
from  the  fact  that  the  spatial  model  for  the  spread  signals  is  stochastic  and  requires  a 
sufficient  number  of  snapshots  to  build  up  the  appropriate  statistics  in  the  data  as  was 
discussed  in  Section  6.2.  For  a  single  point-source  signal,  a  single  snapshot  is  sufficient, 
but  for  a  spread-source,  the  number  of  snapshots  required  to  achieve  a  given  accuracy 
rises  as  the  spread  region  increases  in  size  (e.g.,  compare  the  accuracies  of  the  two  spread 
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signals  in  Figure  9(b)).  Conversely,  for  a  given  number  of  snapshots,  accuracy  degrades 
as  the  spread  region  increases. 

The  effect  of  the  model  uncertainty  can  be  seen  in  the  accuracy  results  for  the  SNR 
=  30  dB  setting.  Compared  to  the  accuracy  results  for  the  SNR  =  15  dB  setting,  the 
improvement  at  the  higher  SNR  is  only  marginal  (e.g.,  compare  this  to  the  improvement 
that  occurs  when  the  SNR  is  increased  from  0  to  15  dB).  Hence  the  main  source  of  error 
at  the  higher  SNR  is  model  uncertainty,  even  for  spreading  as  narrow  as  1°.  For  example, 
from  Figure  14  the  accuracy  for  A ^  =  1°  is  approximately  0.02°  RMS  while  the  accuracy 
for  a  point-source  signal  under  the  same  conditions  is  approximately  0.003°  RMS  (see 
Figure  13). 

The  second  factor,  the  filter  effect,  arises  from  the  fact  that  many  advanced  DF  es¬ 
timators,  such  as  MUSIC  and  SML,  can  be  interpreted  as  techniques  which  use  spatial 
filters  to  suppress  the  signal  content  of  the  input  data  in  the  estimation  procedure  (con¬ 
sider  the  discussion  of  the  whitening  filter  in  Section  3.2).  For  point-source  signals,  this 
requires  creating  a  filter  with  M  notches  corresponding  to  the  M  signals.  These  notches 
are  then  adjusted  so  they  are  located  in  the  same  direction  as  the  signals  themselves. 

The  adjustment  of  the  notch  can  be  made  by  measuring  the  filtered  power  and  then 
adjusting  each  notch  to  minimize  this  power.  Ideally  the  beamwidth  of  each  notch  would 
be  infinitely  narrow  so  that  the  filtered  power  would  be  reduced  only  when  the  notch  was 
placed  exactly  in  the  direction  of  a  signal.  In  reality,  the  width  of  the  notch  is  finite  so 
that  some  noise  also  passes  through  the  filter  which,  in  turn,  affects  the  accuracy  of  the 
adjustment  procedure.  The  width  of  the  notches  is  related  to  the  inverse  of  the  array 
aperture,  hence  an  array  with  a  smaller  aperture  provides  better  accuracy. 

For  spread  signals,  the  beamwidths  of  the  generated  filters  become  much  wider  be¬ 
cause  they  are  designed  to  match  the  spread  widths  of  the  signals.  The  wider  beamwidth 
naturally  results  in  a  greater  amount  of  noise  passing  through  the  filter  with  a  corre¬ 
spondingly  greater  degradation  in  accuracy.  Given  that  the  results  for  SNR  =  30  dB  are 
dominated  by  the  model  uncertainty,  the  results  for  SNR  =  0  dB  are  certainly  dominated 
by  the  filter  effect.  In  this  case,  there  is  little  degradation  in  accuracy  for  spreading  angles 
less  than  4  —  5°  (one  half  the  beamwidth  of  the  array)  compared  to  the  accuracy  of  a 
point-source  signal  (0.1°  RMS).  This  is  due  to  the  limitations  on  the  narrowness  of  the 
filter  beamwidths  discussed  above.  For  greater  spreading  values,  the  change  in  accuracy 
due  to  noise  as  a  function  of  spreading  angle  follows  approximately  the  same  trend  as  for 
the  change  in  accuracy  due  to  model  uncertainty. 
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7.3  Detection  of  a  Weaker  Signal 


As  mentioned  in  the  introduction,  it  is  important  that  a  DF  algorithm  be  able  to 
detect  a  weaker  point-source  signal  in  the  presence  of  stronger  spread-source  signals.  To 
evaluate  the  performance  of  the  SML  algorithm,  simulations  were  run  involving  a  single 
spread-source  with  a  fixed  bearing  and  a  single  point-source  whose  bearing  was  adjusted 
incrementally,  beginning  with  a  large  initial  angular  difference,  until  the  two  bearings 
coincided.  After  each  increment,  the  signal  power  of  the  point-source  was  increased  from 
a  low  value  in  0.5  dB  intervals  until  the  failure  rate  dropped  below  5%  (5  out  of  100  trials). 
The  corresponding  SNR  of  the  point-source  signal  is  defined  here  as  the  threshold  SNR 
and  provides  a  good  indication  of  the  limits  of  detectability  of  the  point-source  signal 
for  the  given  signal  environment.  The  relevant  signal  and  noise  parameters  are  shown  in 
Table  4. 


Table  4:  Signal  Parameters  for  Signal  Detectability  Simulation 


Signal 

<t> 

A* 

Power 

1 

2 

noise 

180° 

adjusted 

30° 

30° 

30° 

0° 

15° 

0° 

0  dB 
adjusted 
-20  dB 

The  results  from  the  simulations  are  shown  in  Figure  15.  In  this  case  the  failure  rate  is 
not  shown  since  it  was  fixed  to  approximately  5%.  From  the  results  it  is  clear  that  when 
the  point-source  and  spread-source  signals  are  well  separated,  25°  or  more  (>  1.2 (/)bw), 
the  presence  of  the  spread  source  signal  has  a  small  but  significant  effect.  For  example, 
the  threshold  SNR  for  a  point-source  without  any  other  signals  present  occurs  at  -12  dB 
(see  Figure  13).  In  the  presence  of  the  spread-source,  the  threshold  SNR  ranges  from  -10 
to  -7  dB  which  is  a  degradation  of  2  to  5  dB.  As  the  separation  between  signals  is  reduced 
from  25°  to  0°,  the  threshold  SNR  increases  dramatically  until  it  reaches  a  maximum  at 
16  dB.  However,  even  with  both  signals  arriving  from  coincident  directions,  detection  of 
both  signals  is  not  only  possible,  but  the  minimum  point-source  power  is  still  4  dB  less 
than  the  spread-source  power. 

The  differences  between  the  wide  and  narrow  separation  cases  can  be  understood  in 
terms  using  a  spatial  filter  to  suppress  the  effects  of  the  spread-source  signal.  In  the  wide 
separation  case,  the  filtering  can  be  accomplished  easily,  leaving  only  the  noise  as  the 
main  source  of  error.  In  the  narrow  separation  case,  it  is  more  and  more  difficult  to  filter 
out  the  spread  signal  independently  of  the  point-source  signal  as  the  separation  decreases. 
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Consequently,  the  spread-source  signal  begins  to  act  as  a  strong  noise  background  and 
the  threshold  increases  accordingly. 

The  accuracy  results  for  the  two  signals  are  relatively  constant  for  wide  spacing  and 
relatively  independent  of  the  presence  of  each  other.  For  example,  in  the  single  signal 
case,  the  accuracies  were  measured  to  be  0.15°  for  the  spread-source  and  no  point  source, 
and  0.7°  for  the  point-source  and  no  spread  source  -  worse  for  the  point-source  due  to  the 
much  lower  SNR.  For  narrower  spacings,  the  threshold  signal  power  of  the  point-source 
signal  increases,  degrading  the  accuracy  of  the  spread-source  signal  until  the  separation 
becomes  as  small  as  8°.  Under  8°  separation,  the  accuracy  actually  improves  even  though 
the  signal  power  of  the  point-source  signal  is  still  increasing,  perhaps  because  the  two 
signal  directions  are  too  similar  to  cause  an  appreciable  error.  The  accuracy  of  the  point- 
source  remains  in  the  range  of  0.5°  — 1.5°  for  most  separation  angles  with  the  lowest  errors 
occurring  for  the  smallest  separations. 

The  simulations  were  also  repeated  using  the  MUSIC  algorithm.  The  assumed  number 
of  signal  directions  (a  parameter  also  required  by  this  algorithm)  used  was  six,  since 
occasionally  up  to  five  directions  were  required  to  describe  the  spread-source  leaving  at 
least  one  direction  for  the  point-source.  In  most  cases,  however,  only  three  or  four  signal 
directions  were  required  for  the  spread  signal,  resulting  in  one  or  more  false  direction 
estimates.  This  is  illustrated  in  Figure  16  which  shows  the  results  of  processing  100  trials 
using  MUSIC. 

Generally,  the  necessity  of  using  several  signal  directions  to  describe  a  spread-source 
and  the  problem  of  extraneous  signals,  makes  the  interpretation  of  the  MUSIC  results 
somewhat  problematic.  For  the  sake  of  this  report,  only  the  two  (out  of  six)  estimated 
signal  bearings  closest  to  the  true  signal  bearings  were  used  when  generating  the  statistical 
results.  Additionally,  since  for  separation  angles  less  than  10°  there  was  no  obvious  way 
to  determine  whether  a  peak  was  associated  with  the  spread  signal  or  the  point-source 
signal  (compare  the  examples  shown  in  Figure  17(a)  and  (b)),  no  statistical  results  for 
MUSIC  were  calculated  for  these  angular  separations.  The  failure  and  accuracy  results 
are  shown  as  the  dashed  lines  in  Figure  15. 

Comparing  the  results  for  the  SML  and  MUSIC  algorithms,  for  all  signal  separations 
measured,  SML  performed  better.  In  terms  of  threshold  performance,  the  differences 
were  no  more  than  2.5  dB  for  the  wider  signal  separations,  but  are  as  large  as  7.5  dB  for 
the  narrower  separations.  Additionally,  for  angular  separations  less  than  10°,  the  SML 
algorithm  was  still  able  to  produce  results  which  could  be  unambiguously  associated  with 
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(a) 


(b) 


(c) 


Figure  15:  Effect  of  angular  spacing  on  the  ability  to  detect  a  weaker  point-source  signal 
in  the  presence  of  a  stronger  spread-source  signal  showing  (a)  the  detection  threshold 
SNR  for  the  weaker  point-source  signal,  (b)  accuracy  of  the  spread-source  estimates  at 
threshold,  and  (c)  accuracy  of  the  point-source  estimates  at  threshold.  The  solid  lines 
represent  the  SML  results,  and  the  dashed  lines  represent  the  MUSIC  results. 
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Figure  16:  MUSIC  azimuth  bearing  estimates  for  100  trials  assuming  6  signal  directions. 
Bearing  estimates  of  the  spread-source  signal  at  180°  and  the  point  source  signal  at  90° 
are  clearly  evident  along  with  various  false  signal  bearing  estimates  (particularly  at  222°). 

the  two  signals  using  either  signal  power  or  spread  estimates,  whereas  with  the  MUSIC 
algorithm  this  was  not  possible  as  previously  discussed. 

The  accuracy  in  estimating  the  point-source  signal  direction  was  the  same  for  both 
algorithms.  However,  since  the  SML  algorithm  has  a  lower  threshold,  this  implies  that 
tested  at  the  same  SNR  (e.g.,  the  MUSIC  threshold)  the  SML  algorithm  is  more  accurate 
than  MUSIC.  The  accuracy  of  MUSIC  for  estimating  the  direction  of  the  spread-source 
signal  was  very  poor,  highlighting  the  difficulty  of  estimating  spread-source  signals  using 
a  point-source  model. 

Generally  the  results  show  that  for  signal  environments  with  spread-source  signals, 
improved  modeling  leads  to  significantly  better  performance. 


7.4  Performance  using  High  Latitude  Off-Air  Data 

The  ultimate  test  for  the  SML  algorithm  is  against  the  data  for  which  it  was  originally 
developed,  namely  actual  high  latitude  HF  measurement  data.  To  this  end,  an  off-air 
measurement  data  set  was  processed  which  had  been  collected  on  January  23,  1996,  using 
the  Vortex  DF  system  at  CFB  Alert  which  is  located  on  the  northern  tip  of  Ellesmere 
Island  in  Northern  Canada  (82.50°  N,  62.35°  W).  The  received  signal  was  from  the  CZD 
transmitter  located  further  south  at  Iqaluit  (63.45°  N,  68.30°  W)  and  transmitting  at  a 
frequency  of  9.292  Mhz.  The  great  circle  signal  bearing  (azimuth)  was  188.5°  measured 
clockwise  from  north.  The  signal  itself  was  composed  of  a  15  second  tone  followed  by  a 
15  second  Morse  code  call  sign  repeated  continuously  for  a  total  period  of  25  minutes. 
This  was,  in  turn,  followed  by  a  5  minute  sounder  slot  during  which  time  there  was  no 


(a) 


(b) 

Figure  17:  Detection  of  a  weaker  point-source  signal  in  the  presence  of  a  stronger  spread- 
source  signal  using  MUSIC.  The  azimuth  spectrum  is  shown  for  an  elevation  angle  of 
*l>  =  30°.  The  simulation  parameters  for  the  noise  and  spread-source  are  listed  in  Table  4, 
and  the  point-source  parameters  were  (a)  <j>  =  150°  with  a  signal  power  of -7.5  dB  and  (b) 
0  =  172°  with  a  signal  power  of  14.5  dB.  When  the  separation  between  the  spread-source 
and  point-source  is  too  narrow,  as  in  (b),  there  is  confusion  as  to  which  spectral  peak 
belongs  to  which  signal. 
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transmission.  This  whole  sequence  was  then  repeated  over  and  over  again. 

The  Vortex  system  consisted  of  an  array  of  12  antennas  connected  to  a  12  channel 
receiver.  The  antenna  array  arrangement,  which  is  as  shown  in  Figure  18,  utilized  8 
elevated  feed  monopole  antennas  from  the  inner  ring  of  a  Pusher  array  (a  circular  array 
with  24  antennas)  and  4  outlying  antennas  (also  elevated  feed  monopole  antennas)  used 
to  increase  the  array  aperture. 


Figure  18:  Three  dimensional  view  of  the  Vortex  antenna  array.  Each  grid  square  has  a 
dimension  of  1A  x  1A. 

The  Vortex  receivers  were  used  to  downconvert  the  input  signals  from  HF  to  2.5  kHz 
with  a  filtered  bandwidth  of  3.5  kHz.  The  downconverted  signals  were  then  digitized  at 
a  rate  of  10  kHz. 

To  generate  the  covariance  estimates,  an  FFT  was  first  performed  on  each  block  of 
12  x  16000  data  points  and  only  the  positive  frequency  data  corresponding  to  3.8  to  4.2 
kHz  were  retained.  This  served  the  dual  purpose  of  converting  the  data  to  IQ  format 
and  suppressing  interference  due  to  noise  and  other  unintended  HF  signals.  The  data 
covariance  matrix  was  then  formed  directly  from  this  frequency  domain  data. 

The  particular  time  of  the  data  collection  was  a  very  active  period  in  terms  of  iono¬ 
spheric  disturbances.  Previous  analysis  using  the  MUSIC  algorithm  [9]  found  that  there 
were  very  large  bearing  swings  over  time.  A  high  degree  of  scatter  in  the  measured  bear¬ 
ings  was  also  observed,  which  is  indicative  of  spread  signals.  The  results  of  reprocessing 
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the  data  as  described  here  and  using  the  MUSIC  algorithm  are  shown  in  Figure  19.  For 
this  processing,  the  number  of  signals  parameter  was  set  to  7  which  was  very  likely  an 
overestimate  of  the  number  of  signal  directions  required,  but  which  was  far  less  detri¬ 
mental  to  the  overall  accuracy  than  underestimating  this  number  would  have  been.  For 
display  purposes,  only  the  bearings  corresponding  to  the  four  largest  peaks  in  the  MUSIC 
spectrum  are  plotted  in  the  figure.  Additionally,  thresholding  was  applied  to  reject  peaks 
in  the  MUSIC  spectrum  which  were  clearly  due  to  noise. 
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Figure  19:  MUSIC  results  for  high  latitude  off-air  data  collected  on  January  23rd  1996. 
Only  the  directions  corresponding  to  the  four  highest  peaks  in  the  MUSIC  spectrum  are 
shown  for  (a)  azimuth,  and  (b)  elevation. 

Even  limiting  the  number  of  signal  directions  displayed  leads  to  some  confusion  in  the 
interpretation  of  the  results,  particularly  in  elevation.  To  provide  some  additional  insight, 
the  results  have  been  replotted  in  Figure  20,  except  this  time  showing  only  the  direction 
of  the  signal  with  the  largest  peak.  It  is  quite  evident  from  these  results  that  the  azimuth 


73 


bearing  was  increasing  over  time  and  that  there  was  considerable  scatter  in  the  results. 


(a) 


(b) 


Figure  20:  MUSIC  results  for  high  latitude  off-air  data  collected  on  January  23rd  1996. 
Only  the  direction  corresponding  to  the  largest  peak  in  the  MUSIC  spectrum  is  shown 
for  (a)  azimuth,  and  (b)  elevation. 

Before  discussing  the  SML  results,  a  comment  is  in  order  about  the  block  size.  Al¬ 
though  a  block  size  of  16000  (representing  1.6  seconds)  would  appear  to  be  sufficient  for 
SML  processing  purposes,  this  does  not  take  into  consideration  the  requirement  for  sample 
to  sample  decorrelation.  The  time  required  to  achieve  this  decorrelation  is  related  to  the 
Doppler  spreading  of  the  signal,  and  can  be  approximated  by  r  =  1  /{Doppler  Spread). 
Values  of  about  20  Hz  were  observed  for  the  collection  period  discussed  here  [9]  giving  a 
value  of  r  =  50  mS.  Hence  for  1.6  seconds  of  data,  this  corresponds  to  an  effective  block 
size  of  K  —  32.  Although  a  higher  value  would  be  better  from  the  accuracy  standpoint, 
this  value  was  not  considered  unreasonable. 

In  processing  the  actual  data,  the  estimated  signal  azimuth  and  elevation  spread  values 
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(a) 


(b) 


Figure  21:  SML  estimated  bearing  results  for  high  latitude  off-air  data  collected  on 
January  23rd  1996.  The  results  are  shown  for  (a)  azimuth,  and  (b)  elevation. 
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were  normally  found  to  be  less  than  60°  except  for  signals  with  low  SNR.  Given  the 
likelihood  that  the  very  large  spread  estimates  (>  60°)  were  errors  due  to  noise,  the 
spread  values  used  by  the  SML  algorithm  were  limited  to  the  range  from  0°  to  60°  to 
avoid  unnecessarily  increasing  the  computational  requirements. 

The  results  of  processing  the  data  using  the  SML  algorithm  are  shown  in  Figure  21. 
Since  the  number  of  signals  was  unknown  in  this  case,  the  data  was  processed  assuming 
M  =  3  signals  while  also  retaining  the  results  for  the  intermediate  cases  of  M  =  0,  M  =  1, 
and  M  =  2.  Bearing  estimates  for  data  block  k  were  then  compared  to  data  blocks  k  —  1 
and  k  +  1  for  the  M  =  3  case.  Real  signals  were  counted  as  those  whose  estimates  varied 
by  less  than  one  beamwidth  over  the  span  of  the  three  blocks.  Additionally,  estimated 
signal  bearings  were  only  counted  if  their  corresponding  estimated  SNR  exceeded  -5  dB. 
Estimates  not  meeting  both  criteria  were  assumed  to  be  due  to  noise  or  other  error 
mechanisms.  The  results  for  a  given  data  block  were  then  taken  according  to  the  case 
M  —  m  where  m  was  the  number  of  real  signals  counted. 

The  results  for  the  SML  algorithm  are  considerably  more  informative  than  those  of  the 
MUSIC  algorithm,  with  a  number  of  different  “tracks”  evident.  The  track  corresponding 
to  the  signal  with  most  power  (dominant)  is  shown  in  Figure  22  which  is  similar  to  the 
results  in  Figure  20  except  that  there  is  less  scatter  in  the  estimates.  The  estimated  signal 
spread  for  the  dominant  signal  is  shown  in  Figure  23.  Interestingly,  given  the  observed 
trend  in  the  azimuth  spread,  the  lowest  spread  occurs  during  the  period  8:48  to  8:52  UT 
when  the  azimuth  signal  direction  corresponded  to  the  great  circle  bearing. 

The  estimated  power  of  the  dominant  signal  and  of  the  noise  are  shown  in  Figure  24. 
The  drop-outs  in  the  signal  power  results  (the  measurements  where  signal  power  dropped 
to  approximately  0  dB  after  being  at  higher  levels  in  previous  measurements)  correspond 
to  measurements  when  the  transmitter  was  off  for  all  or  part  of  the  measurement  period. 
Inspecting  the  estimated  noise  power  results  reveals  that  modeling  of  the  data  was  less 
than  ideal  as  the  noise  power  was  correlated  with  the  signal  power  resulting  in  noise 
power  estimates  up  to  10  dB  higher  than  the  true  power  level  (which  could  be  measured 
during  periods  when  the  transmitter  was  off).  The  cause  of  the  modeling  error  is  probably 
due  to  array  calibration  problems.  For  example,  mutual  coupling  will  have  had  an  effect 
on  the  antennas  selected  from  the  inner  circular  ring  of  the  Pusher  array,  but  this  was 
not  accounted  for  in  the  processing.  This  kind  of  miscalibration  would  then  “warp”  the 
apparent  spatial  distribution  of  a  received  spread  signal  causing  errors. 

The  other  signal  tracks  shown  in  Figure  21  were  due  to  signals  with  low  SNR  (typically 


76 


less  than  0  dB)  which  leaves  open  the  possibility  that  some  of  these  tracks  may  have 
been  processing  artifacts  due  to  the  calibration  problems  already  mentioned.  Despite 
these  problems,  however,  the  SML  algorithm  still  produced  results  which  had  much  less 
scatter  than  MUSIC.  Additionally,  the  parameter  estimates  from  the  SML  algorithm  could 
be  more  easily  used  to  determine  the  number  of  signals  by  comparing  sequential  signal 
parameter  estimates,  to  filter  out  wild  bearings  using  SNR,  and  to  provide  a  useful  quality 
estimate  based  on  both  SNR  and  signal  spreading. 


(a) 


(b) 

Figure  22:  SML  results  for  high  latitude  off-air  data  collected  on  January  23rd  1996. 
Only  the  estimated  bearing  corresponding  to  the  dominant  signal  is  shown  for  (a)  azimuth, 
and  (b)  elevation. 
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(a) 


(b) 


Figure  23:  SML  spreading  results  for  the  dominant  signal  showing  the  (a)  estimated 
azimuth  spreading,  and  (b)  the  estimated  elevation  spreading.  Spreading  estimates  were 
limited  to  the  range  from  0°  to  60°. 
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Universal  Time 


(b) 

Figure  24:  SML  results  showing  (a)  the  estimated  power  of  the  dominant  signal  and  (b) 
the  estimated  noise  power.  Both  power  levels  are  shown  relative  to  the  true  noise  power. 


8.0  CONCLUSIONS  AND  RECOMMENDATIONS 


In  this  report  a  new  direction  finding  estimator  based  on  maximum  likelihood  princi¬ 
ples  was  introduced.  It  is  distinguished  from  other  maximum  likelihood  DF  estimators  by 
the  fact  that  it  deals  with  signals  that  are  spatially  spread  in  terms  of  both  the  received 
azimuth  and  elevation  angles.  For  this  reason  it  is  called  the  Spread  Maximum  Likelihood 
(SML)  algorithm. 

The  implementation  of  the  actual  algorithm  relies  loosely  on  an  approach  called  the  al¬ 
ternating  projection  maximization  method  [11].  The  main  differences  are  that  it  performs 
more  checking  of  intermediate  solutions  (for  improved  accuracy)  and  it  uses  a  gradient  as¬ 
cent  technique  as  the  search  engine.  The  gradient  ascent  technique  has  also  been  modified 
from  the  normal  approach  to  improve  convergence  speed. 

Testing  of  the  SML  algorithm  to  provide  guidelines  on  the  various  control  parameters 
was  performed  through  simulation.  These  parameters  included  the  grid  spacing  used  for 
modeling  the  spread  signals,  the  data  block  size  (the  amount  of  data  required  to  generate 
the  data  covariance  matrix),  the  spread  signal  shape,  and  the  assumed  number  of  signals. 
The  testing  showed  that  these  parameters  could  be  varied  for  a  range  of  values  over 
which  the  performance  of  the  SML  algorithm  was  either  almost  unaffected,  or  was  very 
predictable. 

Testing  of  the  SML  algorithm  was  also  carried  out  through  simulation  to  evaluate  its 
performance  as  a  function  of  various  conditions  of  the  signal  environment.  This  included 
the  effect  of  noise,  signal  spreading,  and  the  detection  of  a  weaker  signal  in  the  presence  of 
a  stronger  signal.  In  all  cases,  the  performance  of  the  SML  algorithm  was  predictable,  and 
as  good  as  or  better  (particularly  for  spatially  spread  signals)  than  the  MUSIC  algorithm 
(a  modern  superresolution  DF  algorithm). 

Finally,  testing  was  performed  using  off-air  data  collected  at  CFS  Alert  which  is  located 
in  the  Arctic.  Again  the  results  showed  that  the  SML  algorithm  outperformed  MUSIC 
in  terms  of  the  variance  in  the  bearing  estimates.  The  SML  estimates  also  included  the 
amount  of  signal  spreading  (MUSIC  does  not)  which  was  on  the  order  of  tens  of  degrees  in 
azimuth.  From  an  understanding  of  the  mechanisms  giving  rise  to  signal  spreading,  these 
measurements  alone  indicate  that  the  signals  were  scattered  from  moving  patches.  The 
observed  trend  of  the  bearing  estimates  over  time  essentially  confirmed  this.  Although  no 
evidence  of  sporadic-E  propagation  was  found  for  this  particular  data  set,  it  was  observed 


that  the  lowest  amount  of  signal  spreading  occurred  when  the  estimated  azimuth  bearing 
coincided  with  the  great  circle  bearing. 

The  main  advantage  of  the  SML  algorithm  is  that,  compared  to  conventional  and 
other  modern  superresolution  techniques,  it  better  models  the  high  latitude  HF  signal 
environment.  This  leads  to  more  accurate  bearing  estimates,  less  confusion  in  the  inter¬ 
pretation  of  the  results  (MUSIC  interprets  a  spread  signal  as  a  multitude  of  point-source 
signals),  and  a  greater  ability  to  detect  a  weaker  point-source  signal  in  the  presence  of 
stronger  spatially  spread  signals.  This  last  point  is  important  since  a  point-source  signal 
will  generally  yield  a  bearing  estimate  close  to  the  great  circle  bearing  while  a  spread 
signal  will  not.  Additionally,  the  ability  of  the  SML  algorithm  to  measure  spatial  signal 
spreading  provides  valuable  information  on  the  reliability  of  the  corresponding  bearing 
estimate. 

The  main  disadvantage  of  the  SML  algorithm  is  that  it  is  computationally  very  in¬ 
tensive  which  makes  it  an  order  of  magnitude  or  more  slower  than  other  approaches.  For 
more  benign  propagation  conditions,  current  superresolution  techniques  such  as  MUSIC 
might  be  more  appropriate  for  use. 

There  are  two  main  recommendations  stemming  from  the  work  detailed  in  this  report. 
The  first  recommendation  is  that  for  high  latitude  operational  HF  DF  sites,  given  its 
benefits,  the  SML  algorithm  should  be  employed  in  some  form.  For  example,  given  the 
computational  requirements  of  the  SML  algorithm,  this  algorithm  could  be  used  in  con¬ 
junction  with  faster  approaches  such  as  MUSIC.  The  faster  approach  could  be  used  until 
spread  conditions  are  suspected,  such  as  when  a  number  of  signal  bearings  are  estimated 
from  similar  directions,  at  which  point,  the  SML  algorithm  would  then  be  used.  The 
second  recommendation  is  that  future  research  should  focus  on  developing  realtime  im¬ 
plementations  of  the  SML  algorithm.  Coupled  with  the  ever  increasing  processing  speed 
of  modern  computers,  a  realtime  implementation  could  be  a  practical  reality  within  a  few 
years. 
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